#include <FilterIntegral.h>

Inheritance diagram for Nektar::SolverUtils::FilterIntegral:

Public Member Functions
SOLVER_UTILS_EXPORT	FilterIntegral (const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< EquationSystem > &pEquation, const ParamMap &pParams)
	Constructs the integral filter and parses filter options, opens file. More...

SOLVER_UTILS_EXPORT	~FilterIntegral () override=default
	Default destructor. More...

Public Member Functions inherited from Nektar::SolverUtils::Filter
SOLVER_UTILS_EXPORT	Filter (const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_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 ()

SOLVER_UTILS_EXPORT std::string	SetupOutput (const std::string ext, const ParamMap &pParams)

SOLVER_UTILS_EXPORT std::string	SetupOutput (const std::string ext, const std::string inname)

SOLVER_UTILS_EXPORT void	SetUpdateOnInitialise (bool flag)

Static Public Member Functions
static FilterSharedPtr	create (const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< EquationSystem > &pEquation, const std::map< std::string, std::string > &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) final
	Parse composite list and geometric entities to integrate over. More...

void	v_Update (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) final
	Calculate integral over requested composites. More...

void	v_Finalise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) final
	This is a time-dependent filter. More...

bool	v_IsTimeDependent () final
	Returns true as filter depends on time. More...

Protected Member Functions inherited from Nektar::SolverUtils::Filter
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

virtual SOLVER_UTILS_EXPORT std::string	v_SetupOutput (const std::string ext, const ParamMap &pParams)

virtual SOLVER_UTILS_EXPORT std::string	v_SetupOutput (const std::string ext, const std::string inname)

Private Attributes
size_t	m_index = 0

size_t	m_outputFrequency
	Frequency to write to output data file in timesteps. More...

size_t	m_numVariables
	Number of fields to perform integral on. More...

std::ofstream	m_outFile
	Out file. More...

LibUtilities::CommSharedPtr	m_comm
	Global communicator. More...

std::vector< std::string >	m_splitCompString
	Vector of composite IDs as a single string. More...

std::vector< std::vector< unsigned int > >	m_compVector
	Vector of vector of composites IDs as integers. More...

std::map< size_t, size_t >	m_geomElmtIdToExpId
	Mapping from geometry ID to expansion ID. More...

std::map< int, std::vector< std::pair< LocalRegions::ExpansionSharedPtr, int > > >	m_compExpMap
	Map of composite ID to vector of expansions an face/edge local ID. More...

Friends
class	MemoryManager< FilterIntegral >

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

bool	m_updateOnInitialise = true

Detailed Description

Definition at line 42 of file FilterIntegral.h.

Constructor & Destructor Documentation

◆ FilterIntegral()

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

Constructs the integral filter and parses filter options, opens file.

Initialise the filter and parse the session file parameters.

Definition at line 52 of file FilterIntegral.cpp.

    : Filter(pSession, pEquation)
{
    std::string outName;
    std::string ext = ".int";
    outName         = Filter::SetupOutput(ext, pParams);
 
    // Composites (to calculate integrals on)
    auto it = pParams.find("Composites");
    ASSERTL0(it != pParams.end(), "Missing parameter 'Composites'.");
    ASSERTL0(it->second.length() > 0, "Empty parameter 'Composites'.");
 
    std::vector<std::string> splitComposite;
    boost::split(m_splitCompString, it->second, boost::is_any_of(","),
                 boost::token_compress_on);
 
    for (auto &comp : m_splitCompString)
    {
        boost::trim(comp);
        size_t first   = comp.find_first_of('[') + 1;
        size_t last    = comp.find_last_of(']') - 1;
        auto tmpString = comp.substr(first, last - first + 1);
 
        std::vector<unsigned int> tmpVec;
        bool parseGood = ParseUtils::GenerateSeqVector(tmpString, tmpVec);
 
        ASSERTL0(parseGood && !tmpVec.empty(),
                 "Unable to read composite regions index range for "
                 "FilterIntegral: " +
                     comp);
 
        m_compVector.emplace_back(tmpVec);
    }
 
    // OutputPrecision
    size_t precision;
    it = pParams.find("OutputPrecision");
    if (it == pParams.end())
    {
        precision = 7;
    }
    else
    {
        ASSERTL0(it->second.length() > 0, "Empty parameter 'OutputPrecision'.");
        precision = std::stoi(it->second);
    }
 
    // Lock equation system pointer
    auto equationSys = m_equ.lock();
    ASSERTL0(equationSys, "Weak pointer expired");
 
    m_numVariables = equationSys->GetNvariables();
 
    m_comm = pSession->GetComm();
    if (m_comm->GetRank() == 0)
    {
        m_outFile.open(outName);
        ASSERTL0(m_outFile.good(), "Unable to open: '" + outName + "'");
        m_outFile.setf(std::ios::scientific, std::ios::floatfield);
        m_outFile.precision(precision);
        m_outFile << "#Time";
 
        for (auto &compName : m_splitCompString)
        {
            for (size_t j = 0; j < m_numVariables; ++j)
            {
                std::string varName = equationSys->GetVariable(j);
                m_outFile << " " + compName + "_" + varName + "_integral";
            }
        }
        m_outFile << std::endl;
    }
 
    // OutputFrequency
    it = pParams.find("OutputFrequency");
    if (it == pParams.end())
    {
        m_outputFrequency = 1;
    }
    else
    {
        ASSERTL0(it->second.length() > 0, "Empty parameter 'OutputFrequency'.");
        LibUtilities::Equation equ(m_session->GetInterpreter(), it->second);
        m_outputFrequency = round(equ.Evaluate());
    }
}

References ASSERTL0, Nektar::LibUtilities::Equation::Evaluate(), Nektar::ParseUtils::GenerateSeqVector(), m_comm, m_compVector, Nektar::SolverUtils::Filter::m_equ, m_numVariables, m_outFile, m_outputFrequency, Nektar::SolverUtils::Filter::m_session, m_splitCompString, and Nektar::SolverUtils::Filter::SetupOutput().

◆ ~FilterIntegral()

SOLVER_UTILS_EXPORT Nektar::SolverUtils::FilterIntegral::~FilterIntegral ( )

overridedefault

Default destructor.

Member Function Documentation

◆ create()

static FilterSharedPtr Nektar::SolverUtils::FilterIntegral::create	(	const LibUtilities::SessionReaderSharedPtr &	pSession,
		const std::shared_ptr< EquationSystem > &	pEquation,
		const std::map< std::string, std::string > &	pParams
	)

inlinestatic

Creates an instance of this class.

Definition at line 48 of file FilterIntegral.h.

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

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

◆ v_Finalise()

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

finalprotectedvirtual

This is a time-dependent filter.

Closes the output data file

Parameters

pFields	Field data
time	Current time

Implements Nektar::SolverUtils::Filter.

Definition at line 379 of file FilterIntegral.cpp.

{
    if (m_comm->GetRank() == 0)
    {
        m_outFile.close();
    }
}

References m_comm, and m_outFile.

◆ v_Initialise()

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

finalprotectedvirtual

Parse composite list and geometric entities to integrate over.

Initialises the integral filter and stores the composite expansions in m_compExpMap for all composites specified in the XML file

Parameters

pFields	Field data
time	Current time

Implements Nektar::SolverUtils::Filter.

Definition at line 144 of file FilterIntegral.cpp.

{
 
    // Create map from element ID -> expansion list ID
    auto expList   = pFields[0]->GetExp();
    auto meshGraph = pFields[0]->GetGraph();
 
    for (size_t i = 0; i < expList->size(); ++i)
    {
        auto exp                                           = (*expList)[i];
        m_geomElmtIdToExpId[exp->GetGeom()->GetGlobalID()] = i;
    }
 
    // Create a map from geom ID -> trace expansion list ID
    std::map<size_t, size_t> geomIdToTraceId;
    auto trace = pFields[0]->GetTrace()->GetExp();
    for (size_t i = 0; i < trace->size(); ++i)
    {
        auto exp                                       = (*trace)[i];
        geomIdToTraceId[exp->GetGeom()->GetGlobalID()] = i;
    }
 
    // Get comp list dimension from first composite & element
    auto composites = pFields[0]->GetGraph()->GetComposites();
    size_t meshDim  = pFields[0]->GetGraph()->GetMeshDimension();
 
    for (int i = 0; i < m_compVector.size(); ++i)
    {
        // Check composite is present in the rank
        if (composites.find(m_compVector[i][0]) == composites.end())
        {
            continue;
        }
 
        std::vector<std::shared_ptr<SpatialDomains::Geometry>> geomVec =
            composites[m_compVector[i][0]]->m_geomVec;
        size_t dim =
            composites[m_compVector[i][0]]->m_geomVec[0]->GetShapeDim();
 
        // Vector of all geometry IDs contained within the composite list
        std::vector<size_t> compGeomIds;
        for (auto compNum : m_compVector[i])
        {
            ASSERTL0(composites.find(compNum) != composites.end(),
                     "In FilterIntegral defined composite C[" +
                         std::to_string(compNum) +
                         "] does not exist in the mesh.")
 
            auto compGeom = composites[compNum]->m_geomVec;
 
            for (auto &geom : compGeom)
            {
                compGeomIds.emplace_back(geom->GetGlobalID());
            }
 
            // Only check first element in each comp for dimension
            ASSERTL0(
                dim == compGeom[0]->GetShapeDim(),
                "Differing geometry dimensions specified in FilterIntegral '" +
                    m_splitCompString[i] + "'.");
        }
 
        std::vector<std::pair<LocalRegions::ExpansionSharedPtr, int>>
            tmpCompExp(compGeomIds.size());
 
        // If dimension of composite == dimension of mesh then we only need the
        // expansion of the element
        if (dim == pFields[0]->GetShapeDimension())
        {
            for (size_t j = 0; j < compGeomIds.size(); ++j)
            {
                tmpCompExp[j] = std::make_pair(
                    (*expList)[m_geomElmtIdToExpId[compGeomIds[j]]], -1);
            }
        }
        // however if the dimension is less we need the expansion of the element
        // containing the global composite geometry and the face/edge local ID
        // within that. 3D mesh -> 2D, 2D -> 1D.
        // @TODO: Restructure with the new dimension independent functions
        //        and check all is correct with this filter.
        else if (meshDim == 3 && dim == 2)
        {
            for (size_t j = 0; j < compGeomIds.size(); ++j)
            {
                LocalRegions::ExpansionSharedPtr exp =
                    (*trace)[geomIdToTraceId[compGeomIds[j]]];
                LocalRegions::Expansion2DSharedPtr exp2D =
                    std::dynamic_pointer_cast<LocalRegions::Expansion2D>(exp);
 
                LocalRegions::ExpansionSharedPtr leftAdjElmtExp =
                    std::dynamic_pointer_cast<LocalRegions::Expansion>(
                        exp2D->GetLeftAdjacentElementExp());
                int leftAdjElmtFace = exp2D->GetLeftAdjacentElementTrace();
 
                tmpCompExp[j] = std::make_pair(leftAdjElmtExp, leftAdjElmtFace);
            }
        }
        else if (meshDim == 2 && dim == 1)
        {
            for (size_t j = 0; j < compGeomIds.size(); ++j)
            {
                LocalRegions::ExpansionSharedPtr exp =
                    (*trace)[geomIdToTraceId[compGeomIds[j]]];
                LocalRegions::Expansion1DSharedPtr exp1D =
                    std::dynamic_pointer_cast<LocalRegions::Expansion1D>(exp);
 
                LocalRegions::ExpansionSharedPtr leftAdjElmtExp =
                    std::dynamic_pointer_cast<LocalRegions::Expansion>(
                        exp1D->GetLeftAdjacentElementExp());
                int leftAdjElmtEdge = exp1D->GetLeftAdjacentElementTrace();
 
                tmpCompExp[j] = std::make_pair(leftAdjElmtExp, leftAdjElmtEdge);
            }
        }
        else
        {
            ASSERTL0(false,
                     "FilterIntegral: Only composite dimensions equal to or "
                     "one lower than the mesh dimension are supported.")
        }
 
        m_compExpMap[i] = tmpCompExp;
    }
    if (m_updateOnInitialise)
    {
        v_Update(pFields, time);
    }
}

References ASSERTL0, m_compExpMap, m_compVector, m_geomElmtIdToExpId, m_splitCompString, Nektar::SolverUtils::Filter::m_updateOnInitialise, and v_Update().

◆ v_IsTimeDependent()

bool Nektar::SolverUtils::FilterIntegral::v_IsTimeDependent ( )

finalprotectedvirtual

Returns true as filter depends on time.

Implements Nektar::SolverUtils::Filter.

Definition at line 390 of file FilterIntegral.cpp.

{
    return true;
}

◆ v_Update()

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

finalprotectedvirtual

Calculate integral over requested composites.

Performs the integration on the stored composite expansions and outputs in to the output data file

Parameters

pFields	Field data
time	Current time

Implements Nektar::SolverUtils::Filter.

Definition at line 278 of file FilterIntegral.cpp.

{
    if (m_index++ % m_outputFrequency > 0)
    {
        return;
    }
 
    if (m_comm->GetRank() == 0)
    {
        m_outFile << time;
    }
 
    for (size_t j = 0; j < m_compVector.size(); ++j)
    {
        for (size_t i = 0; i < m_numVariables; ++i)
        {
            Array<OneD, const NekDouble> phys;
            phys = pFields[i]->GetPhys();
 
            NekDouble sum = 0.0;
            NekDouble c   = 0.0;
 
            // Check if composite is on the rank
            if (m_compExpMap.find(j) != m_compExpMap.end())
            {
                auto compExp   = m_compExpMap[j];
                size_t dim     = compExp[0].first->GetGeom()->GetShapeDim();
                size_t meshDim = pFields[i]->GetGraph()->GetMeshDimension();
 
                // Evaluate integral using improved Kahan–Babuška summation
                // algorithm to reduce numerical error from adding floating
                // points
                for (auto &expPair : compExp)
                {
                    NekDouble input = 0;
                    auto exp        = expPair.first;
 
                    if (meshDim == dim)
                    {
                        size_t offset = pFields[i]->GetPhys_Offset(
                            m_geomElmtIdToExpId[exp->GetGeom()->GetGlobalID()]);
                        input = exp->Integral(phys + offset);
                    }
                    else if (meshDim == 3 && dim == 2)
                    {
                        Array<OneD, NekDouble> facePhys;
                        exp->GetTracePhysVals(expPair.second, exp, phys,
                                              facePhys);
                        input =
                            pFields[i]
                                ->GetTrace()
                                ->GetExp(exp->GetGeom()->GetTid(expPair.second))
                                ->Integral(facePhys);
                    }
                    else if (meshDim == 2 && dim == 1)
                    {
                        Array<OneD, NekDouble> edgePhys;
                        exp->GetTracePhysVals(expPair.second, exp, phys,
                                              edgePhys);
                        input =
                            pFields[i]
                                ->GetTrace()
                                ->GetExp(exp->GetGeom()->GetTid(expPair.second))
                                ->Integral(edgePhys);
                    }
                    else
                    {
                        ASSERTL0(false,
                                 "FilterIntegral: Only composite dimensions "
                                 "equal to or one lower than the mesh "
                                 "dimension are supported.")
                    }
 
                    NekDouble t = sum + input;
                    c += fabs(sum) >= fabs(input) ? (sum - t) + input
                                                  : (input - t) + sum;
                    sum = t;
                }
            }
 
            // Sum integral values from all ranks
            Array<OneD, NekDouble> sumArray(1, sum + c);
            m_comm->AllReduce(sumArray, LibUtilities::ReduceSum);
            if (m_comm->GetRank() == 0)
            {
                m_outFile << " " << sumArray[0];
            }
        }
    }
 
    if (m_comm->GetRank() == 0)
    {
        m_outFile << std::endl;
    }
}

References ASSERTL0, m_comm, m_compExpMap, m_compVector, m_geomElmtIdToExpId, m_index, m_numVariables, m_outFile, m_outputFrequency, and Nektar::LibUtilities::ReduceSum.

Referenced by v_Initialise().

Friends And Related Function Documentation

◆ MemoryManager< FilterIntegral >

friend class MemoryManager< FilterIntegral >

friend

Definition at line 1 of file FilterIntegral.h.

Member Data Documentation

◆ className

std::string Nektar::SolverUtils::FilterIntegral::className

static

Initial value:

=
    GetFilterFactory().RegisterCreatorFunction("Integral",
                                               FilterIntegral::create)

Name of the class.

Definition at line 59 of file FilterIntegral.h.

◆ m_comm

LibUtilities::CommSharedPtr Nektar::SolverUtils::FilterIntegral::m_comm

private

Global communicator.

Definition at line 115 of file FilterIntegral.h.

Referenced by FilterIntegral(), v_Finalise(), and v_Update().

◆ m_compExpMap

std::map<int, std::vector<std::pair<LocalRegions::ExpansionSharedPtr, int> > > Nektar::SolverUtils::FilterIntegral::m_compExpMap

private

Map of composite ID to vector of expansions an face/edge local ID.

Definition at line 124 of file FilterIntegral.h.

Referenced by v_Initialise(), and v_Update().

◆ m_compVector

std::vector<std::vector<unsigned int> > Nektar::SolverUtils::FilterIntegral::m_compVector

private

Vector of vector of composites IDs as integers.

Definition at line 119 of file FilterIntegral.h.

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

◆ m_geomElmtIdToExpId

std::map<size_t, size_t> Nektar::SolverUtils::FilterIntegral::m_geomElmtIdToExpId

private

Mapping from geometry ID to expansion ID.

Definition at line 121 of file FilterIntegral.h.

Referenced by v_Initialise(), and v_Update().

◆ m_index

size_t Nektar::SolverUtils::FilterIntegral::m_index = 0

private

Definition at line 107 of file FilterIntegral.h.

Referenced by v_Update().

◆ m_numVariables

size_t Nektar::SolverUtils::FilterIntegral::m_numVariables

private

Number of fields to perform integral on.

Definition at line 111 of file FilterIntegral.h.

Referenced by FilterIntegral(), and v_Update().

◆ m_outFile

std::ofstream Nektar::SolverUtils::FilterIntegral::m_outFile

private

Out file.

Definition at line 113 of file FilterIntegral.h.

Referenced by FilterIntegral(), v_Finalise(), and v_Update().

◆ m_outputFrequency

size_t Nektar::SolverUtils::FilterIntegral::m_outputFrequency

private

Frequency to write to output data file in timesteps.

Definition at line 109 of file FilterIntegral.h.

Referenced by FilterIntegral(), and v_Update().

◆ m_splitCompString

std::vector<std::string> Nektar::SolverUtils::FilterIntegral::m_splitCompString

private

Vector of composite IDs as a single string.

Definition at line 117 of file FilterIntegral.h.

Referenced by FilterIntegral(), and v_Initialise().

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Private Attributes

Friends

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ FilterIntegral()

◆ ~FilterIntegral()

Member Function Documentation

◆ create()

◆ v_Finalise()

◆ v_Initialise()

◆ v_IsTimeDependent()

◆ v_Update()

Friends And Related Function Documentation

◆ MemoryManager< FilterIntegral >

Member Data Documentation

◆ className

◆ m_comm

◆ m_compExpMap

◆ m_compVector

◆ m_geomElmtIdToExpId

◆ m_index

◆ m_numVariables

◆ m_outFile

◆ m_outputFrequency

◆ m_splitCompString