#include <PreconCfsBRJ.h>

Inheritance diagram for Nektar::PreconCfsBRJ:

Public Member Functions
	PreconCfsBRJ (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &vComm)

	~PreconCfsBRJ () override=default

Public Member Functions inherited from Nektar::PreconCfs
	PreconCfs (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &vComm)

virtual	~PreconCfs ()=default

void	InitObject ()

void	DoPreconCfs (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &flag)

void	BuildPreconCfs (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, const Array< OneD, NekDouble > > &intmp, const NekDouble time, const NekDouble lambda)

bool	UpdatePreconMatCheck (const Array< OneD, const NekDouble > &res, const NekDouble dtLambda)

void	SetOperators (const NekPreconCfsOperators &in)

Static Public Member Functions
static PreconCfsSharedPtr	create (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &vComm)
	Creates an instance of this class. More...

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

Protected Member Functions
void	v_InitObject () override

void	v_DoPreconCfs (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &flag) override

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

bool	v_UpdatePreconMatCheck (const Array< OneD, const NekDouble > &res, const NekDouble dtLambda) override

virtual void	v_InitObject ()=0

virtual void	v_DoPreconCfs (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &flag)=0

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

virtual bool	v_UpdatePreconMatCheck (const Array< OneD, const NekDouble > &res, const NekDouble dtLambda)=0

Protected Attributes
int	m_PreconItsStep

int	m_BRJRelaxParam

Array< OneD, Array< OneD, SNekBlkMatSharedPtr > >	m_PreconMatVarsSingle

unsigned int	m_max_nblocks

unsigned int	m_max_nElmtDof

std::vector< simd< NekSingle >, tinysimd::allocator< simd< NekSingle > > >	m_sBlkDiagMat

std::vector< int >	m_inputIdx

Array< OneD, SNekBlkMatSharedPtr >	m_TraceJacSingle

TensorOfArray4D< NekSingle >	m_TraceJacArraySingle

Array< OneD, SNekBlkMatSharedPtr >	m_TraceJacDerivSingle

TensorOfArray4D< NekSingle >	m_TraceJacDerivArraySingle

Array< OneD, Array< OneD, NekSingle > >	m_TraceJacDerivSignSingle

TensorOfArray5D< NekSingle >	m_TraceIPSymJacArraySingle

Protected Attributes inherited from Nektar::PreconCfs
LibUtilities::CommSharedPtr	m_Comm

bool	m_verbose

int	m_spacedim

NekPreconCfsOperators	m_operator

int	m_PreconMatFreezNumb

int	m_PreconTimesCounter

NekDouble	m_DtLambdaPreconMat = -1.0

bool	m_CalcPreconMatFlag = false

Private Member Functions
void	DoNullPrecon (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &flag)

void	PreconBlkDiag (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)

template<typename DataType >
void	MinusOffDiag2Rhs (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const size_t nvariables, const size_t nCoeffs, const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, TensorOfArray3D< NekDouble > &wspTrace, Array< OneD, Array< OneD, DataType > > &wspTraceDataType, const TensorOfArray4D< DataType > &TraceJacArray)

template<typename TypeNekBlkMatSharedPtr >
void	AllocatePreconBlkDiagCoeff (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, Array< OneD, TypeNekBlkMatSharedPtr > > &gmtxarray, const int &nscale=1)

void	AllocateSIMDPreconBlkMatDiag (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
	This function creates the matrix structure for the block diagonal operator. It organizes the way that the matrix. More...

void	AllocateNekBlkMatDig (SNekBlkMatSharedPtr &mat, const Array< OneD, unsigned int > nrow, const Array< OneD, unsigned int > ncol)

Friends
class	MemoryManager< PreconCfsBRJ >

Detailed Description

Block Relaxed(weighted) Jacobi iterative (BRJ) Preconditioner for CFS

Solves a linear system using iterative methods.

Definition at line 49 of file PreconCfsBRJ.h.

Constructor & Destructor Documentation

◆ PreconCfsBRJ()

Nektar::PreconCfsBRJ::PreconCfsBRJ	(	const Array< OneD, MultiRegions::ExpListSharedPtr > &	pFields,
		const LibUtilities::SessionReaderSharedPtr &	pSession,
		const LibUtilities::CommSharedPtr &	vComm
	)

Definition at line 53 of file PreconCfsBRJ.cpp.

    : PreconCfs(pFields, pSession, vComm)
{
    pSession->LoadParameter("PreconItsStep", m_PreconItsStep, 7);
    pSession->LoadParameter("BRJRelaxParam", m_BRJRelaxParam, 1.0);
 
    size_t nvariables = pFields.size();
 
    m_PreconMatVarsSingle =
        Array<OneD, Array<OneD, SNekBlkMatSharedPtr>>(nvariables);
    for (size_t i = 0; i < nvariables; i++)
    {
        m_PreconMatVarsSingle[i] = Array<OneD, SNekBlkMatSharedPtr>(nvariables);
    }
    AllocatePreconBlkDiagCoeff(pFields, m_PreconMatVarsSingle);
 
    AllocateSIMDPreconBlkMatDiag(pFields);
}

References AllocatePreconBlkDiagCoeff(), AllocateSIMDPreconBlkMatDiag(), m_BRJRelaxParam, m_PreconItsStep, and m_PreconMatVarsSingle.

◆ ~PreconCfsBRJ()

Nektar::PreconCfsBRJ::~PreconCfsBRJ ( )

overridedefault

Member Function Documentation

◆ AllocateNekBlkMatDig()

void Nektar::PreconCfsBRJ::AllocateNekBlkMatDig	(	SNekBlkMatSharedPtr &	mat,
		const Array< OneD, unsigned int >	nrow,
		const Array< OneD, unsigned int >	ncol
	)

inlineprivate

Definition at line 236 of file PreconCfsBRJ.h.

    {
        mat =
            MemoryManager<SNekBlkMat>::AllocateSharedPtr(nrow, ncol, eDIAGONAL);
        SNekMatSharedPtr loc_matNvar;
        for (size_t nelm = 0; nelm < nrow.size(); ++nelm)
        {
            int nrowsVars = nrow[nelm];
            int ncolsVars = ncol[nelm];
 
            loc_matNvar = MemoryManager<SNekMat>::AllocateSharedPtr(
                nrowsVars, ncolsVars, 0.0);
            mat->SetBlock(nelm, nelm, loc_matNvar);
        }
    }

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and Nektar::eDIAGONAL.

Referenced by AllocatePreconBlkDiagCoeff(), and v_BuildPreconCfs().

◆ AllocatePreconBlkDiagCoeff()

template<typename TypeNekBlkMatSharedPtr >

void Nektar::PreconCfsBRJ::AllocatePreconBlkDiagCoeff	(	const Array< OneD, MultiRegions::ExpListSharedPtr > &	pFields,
		Array< OneD, Array< OneD, TypeNekBlkMatSharedPtr > > &	gmtxarray,
		const int &	nscale = `1`
	)

private

Definition at line 499 of file PreconCfsBRJ.cpp.

{
 
    size_t nvars  = pFields.size();
    size_t nelmts = pFields[0]->GetNumElmts();
    size_t nelmtcoef;
    Array<OneD, unsigned int> nelmtmatdim(nelmts);
    for (size_t i = 0; i < nelmts; i++)
    {
        nelmtcoef      = pFields[0]->GetExp(i)->GetNcoeffs();
        nelmtmatdim[i] = nelmtcoef * nscale;
    }
 
    for (size_t i = 0; i < nvars; i++)
    {
        for (size_t j = 0; j < nvars; j++)
        {
            AllocateNekBlkMatDig(gmtxarray[i][j], nelmtmatdim, nelmtmatdim);
        }
    }
}

References AllocateNekBlkMatDig().

Referenced by PreconCfsBRJ().

◆ AllocateSIMDPreconBlkMatDiag()

void Nektar::PreconCfsBRJ::AllocateSIMDPreconBlkMatDiag ( const Array< OneD, MultiRegions::ExpListSharedPtr > & pFields )

inlineprivate

This function creates the matrix structure for the block diagonal operator. It organizes the way that the matrix.

by SIMD instructions. The degrees of freedom of each element are reorganized, so they are placed in the correct location to perfom SIMD instructions.

Definition at line 138 of file PreconCfsBRJ.h.

    {
        using vec_t = simd<NekSingle>;
 
        int TotMatLen         = 0;
        int TotLen            = 0;
        const auto nTotElmt   = pFields[0]->GetNumElmts();
        const auto nvariables = pFields.size();
        const auto vecwidth   = vec_t::width;
 
        m_max_nblocks  = 0;
        m_max_nElmtDof = 0;
 
        for (int ne = 0; ne < nTotElmt; ne++)
        {
            const auto nElmtDof = pFields[0]->GetNcoeffs(ne) * nvariables;
            const auto nblocks  = nElmtDof / vecwidth;
            unsigned int totblocks =
                (nElmtDof % vecwidth) ? nblocks + 1 : nblocks;
 
            m_max_nblocks =
                (m_max_nblocks > totblocks) ? m_max_nblocks : totblocks;
            m_max_nElmtDof =
                (m_max_nElmtDof > nElmtDof) ? m_max_nElmtDof : nElmtDof;
 
            TotLen += totblocks * vecwidth;
            TotMatLen += nElmtDof * totblocks;
        }
 
        m_sBlkDiagMat.resize(TotMatLen);
        m_inputIdx.resize(TotLen);
 
        // generate a index list of vector width aware mapping from
        // local coeff storage over all variables to elemental storage
        // over variables
        unsigned int ncoeffs = pFields[0]->GetNcoeffs();
        for (int ne = 0, cnt1 = 0; ne < nTotElmt; ne++)
        {
            const auto nElmtCoeff  = pFields[0]->GetNcoeffs(ne);
            const auto nElmtDof    = nElmtCoeff * nvariables;
            const auto nblocks     = nElmtDof / vecwidth;
            const auto nCoefOffset = pFields[0]->GetCoeff_Offset(ne);
            int i                  = 0;
            int i0                 = 0;
            int inOffset, j;
 
            for (int m = 0; m < nvariables; m++)
            {
                inOffset = m * ncoeffs + nCoefOffset;
 
                if (m && (vecwidth - i0 < nElmtCoeff))
                {
                    // May need to add entries from later variables to
                    // remainder of last variable if the vector width
                    // was not exact multiple of number of elemental
                    // coeffs
                    for (i = 0; i0 < vecwidth; ++i, ++i0)
                    {
                        m_inputIdx[cnt1++] = inOffset + i;
                    }
                }
                else
                {
                    i = 0;
                }
 
                // load up other vectors in variable that fit into vector
                // width
                for (j = 0; (j + 1) * vecwidth < nElmtCoeff - i; ++j)
                {
                    for (i0 = 0; i0 < vecwidth; ++i0)
                    {
                        m_inputIdx[cnt1++] = inOffset + i + j * vecwidth + i0;
                    }
                }
 
                // load up any residual data for this variable
                for (i0 = 0, j = j * vecwidth; j < nElmtCoeff - i; ++j, ++i0)
                {
                    m_inputIdx[cnt1++] = inOffset + i + j;
                }
            }
 
            const auto endwidth = nElmtDof - nblocks * vecwidth;
 
            // fill out rest of index to match vector width with last entry
            if (endwidth)
            {
                for (i0 = endwidth; i0 < vecwidth; ++i0)
                {
                    m_inputIdx[cnt1++] = inOffset + i + j - 1;
                }
            }
            ASSERTL1(cnt1 <= TotLen, "m_inputIdx over extended");
        }
    }

References ASSERTL1, m_inputIdx, m_max_nblocks, m_max_nElmtDof, and m_sBlkDiagMat.

Referenced by PreconCfsBRJ().

◆ create()

static PreconCfsSharedPtr Nektar::PreconCfsBRJ::create	(	const Array< OneD, MultiRegions::ExpListSharedPtr > &	pFields,
		const LibUtilities::SessionReaderSharedPtr &	pSession,
		const LibUtilities::CommSharedPtr &	vComm
	)

inlinestatic

Creates an instance of this class.

Definition at line 55 of file PreconCfsBRJ.h.

    {
        PreconCfsSharedPtr p = MemoryManager<PreconCfsBRJ>::AllocateSharedPtr(
            pFields, pSession, vComm);
        return p;
    }

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

◆ DoNullPrecon()

void Nektar::PreconCfsBRJ::DoNullPrecon	(	const Array< OneD, NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput,
		const bool &	flag
	)

private

Definition at line 288 of file PreconCfsBRJ.cpp.

{
    Vmath::Vcopy(pInput.size(), pInput, 1, pOutput, 1);
}

References Vmath::Vcopy().

Referenced by v_DoPreconCfs().

◆ MinusOffDiag2Rhs()

template<typename DataType >

void Nektar::PreconCfsBRJ::MinusOffDiag2Rhs	(	const Array< OneD, MultiRegions::ExpListSharedPtr > &	pFields,
		const size_t	nvariables,
		const size_t	nCoeffs,
		const Array< OneD, const Array< OneD, NekDouble > > &	inarray,
		Array< OneD, Array< OneD, NekDouble > > &	outarray,
		TensorOfArray3D< NekDouble > &	wspTrace,
		Array< OneD, Array< OneD, DataType > > &	wspTraceDataType,
		const TensorOfArray4D< DataType > &	TraceJacArray
	)

private

Definition at line 377 of file PreconCfsBRJ.cpp.

{
    size_t nTracePts = pFields[0]->GetTrace()->GetNpoints();
    size_t npoints   = pFields[0]->GetNpoints();
    size_t nDim      = m_spacedim;
 
    Array<OneD, Array<OneD, NekDouble>> outpnts(nvariables);
    for (size_t i = 0; i < nvariables; i++)
    {
        outpnts[i] = Array<OneD, NekDouble>(npoints, 0.0);
        pFields[i]->BwdTrans(outarray[i], outpnts[i]);
    }
 
    // Store forwards/backwards space along trace space
    Array<OneD, Array<OneD, NekDouble>> Fwd;
    Array<OneD, Array<OneD, NekDouble>> Bwd;
    Array<OneD, Array<OneD, NekDouble>> FwdFlux;
    Array<OneD, Array<OneD, NekDouble>> BwdFlux;
    TensorOfArray3D<NekDouble> numDerivBwd(nDim);
    TensorOfArray3D<NekDouble> numDerivFwd(nDim);
    size_t indexwspTrace = 0;
    Fwd                  = wspTrace[indexwspTrace], indexwspTrace++;
    Bwd                  = wspTrace[indexwspTrace], indexwspTrace++;
    FwdFlux              = wspTrace[indexwspTrace], indexwspTrace++;
    BwdFlux              = wspTrace[indexwspTrace], indexwspTrace++;
 
    LibUtilities::Timer timer;
    for (size_t i = 0; i < nvariables; ++i)
    {
        timer.Start();
        pFields[i]->GetFwdBwdTracePhys(outpnts[i], Fwd[i], Bwd[i]);
        timer.Stop();
        timer.AccumulateRegion("ExpList::GetFwdBwdTracePhys", 10);
    }
 
    size_t indexwspTraceDataType = 0;
    Array<OneD, Array<OneD, DataType>> Fwdarray(nvariables);
    for (size_t m = 0; m < nvariables; ++m)
    {
        Fwdarray[m] = wspTraceDataType[indexwspTraceDataType],
        indexwspTraceDataType++;
    }
    Array<OneD, DataType> Fwdreslt;
    Fwdreslt = wspTraceDataType[indexwspTraceDataType], indexwspTraceDataType++;
 
    for (size_t m = 0; m < nvariables; ++m)
    {
        for (size_t i = 0; i < nTracePts; ++i)
        {
            Fwdarray[m][i] = DataType(Fwd[m][i]);
        }
    }
    for (size_t m = 0; m < nvariables; ++m)
    {
        Vmath::Zero(nTracePts, Fwdreslt, 1);
        for (size_t n = 0; n < nvariables; ++n)
        {
            for (size_t p = 0; p < nTracePts; ++p)
            {
                Fwdreslt[p] += TraceJacArray[0][m][n][p] * Fwdarray[n][p];
            }
        }
        for (size_t i = 0; i < nTracePts; ++i)
        {
            FwdFlux[m][i] = NekDouble(Fwdreslt[i]);
        }
    }
 
    for (size_t m = 0; m < nvariables; ++m)
    {
        for (size_t i = 0; i < nTracePts; ++i)
        {
            Fwdarray[m][i] = DataType(Bwd[m][i]);
        }
    }
    for (size_t m = 0; m < nvariables; ++m)
    {
        Vmath::Zero(nTracePts, Fwdreslt, 1);
        for (size_t n = 0; n < nvariables; ++n)
        {
            for (size_t p = 0; p < nTracePts; ++p)
            {
                Fwdreslt[p] += TraceJacArray[1][m][n][p] * Fwdarray[n][p];
            }
        }
        for (size_t i = 0; i < nTracePts; ++i)
        {
            BwdFlux[m][i] = NekDouble(Fwdreslt[i]);
        }
    }
 
    for (size_t i = 0; i < nvariables; ++i)
    {
        Vmath::Fill(nCoeffs, 0.0, outarray[i], 1);
        timer.Start();
        pFields[i]->AddTraceIntegralToOffDiag(FwdFlux[i], BwdFlux[i],
                                              outarray[i]);
        timer.Stop();
        timer.AccumulateRegion("ExpList::AddTraceIntegralToOffDiag", 10);
    }
 
    for (size_t i = 0; i < nvariables; ++i)
    {
        for (size_t p = 0; p < nCoeffs; ++p)
        {
            outarray[i][p] =
                -m_DtLambdaPreconMat * outarray[i][p] + inarray[i][p];
        }
    }
}

References Nektar::LibUtilities::Timer::AccumulateRegion(), Vmath::Fill(), Nektar::PreconCfs::m_DtLambdaPreconMat, Nektar::PreconCfs::m_spacedim, CellMLToNektar.cellml_metadata::p, Nektar::LibUtilities::Timer::Start(), Nektar::LibUtilities::Timer::Stop(), and Vmath::Zero().

Referenced by v_DoPreconCfs().

◆ PreconBlkDiag()

void Nektar::PreconCfsBRJ::PreconBlkDiag	(	const Array< OneD, MultiRegions::ExpListSharedPtr > &	pFields,
		const Array< OneD, NekDouble > &	inarray,
		Array< OneD, NekDouble > &	outarray
	)

private

Definition at line 298 of file PreconCfsBRJ.cpp.

{
    unsigned int nvariables = pFields.size();
 
    int nTotElmt = pFields[0]->GetNumElmts();
 
    using vec_t         = simd<NekSingle>;
    const auto vecwidth = vec_t::width;
 
    // vectorized matrix multiply
    std::vector<vec_t, tinysimd::allocator<vec_t>> Sinarray(m_max_nblocks);
    std::vector<vec_t, tinysimd::allocator<vec_t>> Soutarray(m_max_nElmtDof);
 
    alignas(vec_t::alignment) std::array<NekSingle, vec_t::width> tmp;
 
    for (int ne = 0, cnt = 0, icnt = 0, icnt1 = 0; ne < nTotElmt; ne++)
    {
        const auto nElmtCoef = pFields[0]->GetNcoeffs(ne);
        const auto nElmtDof  = nElmtCoef * nvariables;
        const auto nblocks   = (nElmtDof % vecwidth) ? nElmtDof / vecwidth + 1
                                                     : nElmtDof / vecwidth;
 
        // gather data into blocks - could probably be done with a
        // gather call? can be replaced with a gather op when working
        for (int j = 0; j < nblocks; ++j, icnt += vecwidth)
        {
            for (int i = 0; i < vecwidth; ++i)
            {
                tmp[i] = inarray[m_inputIdx[icnt + i]];
            }
 
            Sinarray[j].load(tmp.data());
        }
 
        // Do matrix multiply
        // first block just needs multiplying
        vec_t in = Sinarray[0];
        for (int i = 0; i < nElmtDof; ++i)
        {
            Soutarray[i] = m_sBlkDiagMat[cnt++] * in;
        }
 
        // rest of blocks are multiply add operations;
        for (int n = 1; n < nblocks; ++n)
        {
            in = Sinarray[n];
            for (int i = 0; i < nElmtDof; ++i)
            {
                Soutarray[i].fma(m_sBlkDiagMat[cnt++], in);
            }
        }
 
        // get block aligned index for this expansion
        NekSingle val;
        for (int i = 0; i < nElmtDof; ++i)
        {
            // Get hold of data
            Soutarray[i].store(tmp.data());
 
            // Sum vector width
            val = tmp[0];
            for (int j = 1; j < vecwidth; ++j)
            {
                val += tmp[j];
            }
            // put data into outarray
            outarray[m_inputIdx[icnt1 + i]] = NekDouble(val);
        }
 
        icnt1 += nblocks * vecwidth;
    }
}

References m_inputIdx, m_max_nblocks, m_max_nElmtDof, and m_sBlkDiagMat.

Referenced by v_DoPreconCfs().

◆ v_BuildPreconCfs()

void Nektar::PreconCfsBRJ::v_BuildPreconCfs	(	const Array< OneD, MultiRegions::ExpListSharedPtr > &	pFields,
		const Array< OneD, const Array< OneD, NekDouble > > &	intmp,
		const NekDouble	time,
		const NekDouble	lambda
	)

overrideprotectedvirtual

Copy array into column major blocks of vector width

Implements Nektar::PreconCfs.

Definition at line 188 of file PreconCfsBRJ.cpp.

{
    if (m_PreconItsStep > 0)
    {
        SNekBlkMatSharedPtr PreconMatSingle;
        using vec_t    = simd<NekSingle>;
        int nvariables = pFields.size();
        int nelmts     = pFields[0]->GetNumElmts();
        Array<OneD, unsigned int> matdim(nelmts);
        for (int i = 0; i < nelmts; i++)
        {
            matdim[i] = pFields[0]->GetExp(i)->GetNcoeffs() * nvariables;
        }
        AllocateNekBlkMatDig(PreconMatSingle, matdim, matdim);
 
        m_operator.DoCalcPreconMatBRJCoeff(
            intmp, m_PreconMatVarsSingle, PreconMatSingle, m_TraceJacSingle,
            m_TraceJacDerivSingle, m_TraceJacDerivSignSingle,
            m_TraceJacArraySingle, m_TraceJacDerivArraySingle,
            m_TraceIPSymJacArraySingle);
 
        if (m_verbose && m_Comm->GetRank() == 0)
        {
            cout << "     ## CalcuPreconMat " << endl;
        }
 
        // copy matrix to simd layout
        // load matrix
        int cnt             = 0;
        const auto vecwidth = vec_t::width;
 
        alignas(vec_t::alignment) std::array<NekSingle, vec_t::width> tmp;
 
        for (int ne = 0; ne < nelmts; ne++)
        {
            const auto nElmtDof = matdim[ne];
            const auto nblocks  = nElmtDof / vecwidth;
 
            const NekSingle *mmat =
                PreconMatSingle->GetBlockPtr(ne, ne)->GetRawPtr();
            /// Copy array into column major blocks of vector width
            for (int i1 = 0; i1 < nblocks; ++i1)
            {
                for (int j = 0; j < nElmtDof; ++j)
                {
                    for (int i = 0; i < vecwidth; ++i)
                    {
                        tmp[i] = mmat[j + (i1 * vecwidth + i) * nElmtDof];
                    }
                    // store contiguous vec_t array.
                    m_sBlkDiagMat[cnt++].load(tmp.data());
                }
            }
 
            const auto endwidth = nElmtDof - nblocks * vecwidth;
 
            // end rows that do not fit into vector widths
            if (endwidth)
            {
                for (int j = 0; j < nElmtDof; ++j)
                {
                    for (int i = 0; i < endwidth; ++i)
                    {
                        tmp[i] = mmat[j + (nblocks * vecwidth + i) * nElmtDof];
                    }
 
                    for (int i = endwidth; i < vecwidth; ++i)
                    {
                        tmp[i] = 0.0;
                    }
                    m_sBlkDiagMat[cnt++].load(tmp.data());
                }
            }
        }
    }
 
    m_DtLambdaPreconMat  = lambda;
    m_CalcPreconMatFlag  = false;
    m_PreconTimesCounter = 1;
}

References AllocateNekBlkMatDig(), Nektar::NekPreconCfsOperators::DoCalcPreconMatBRJCoeff(), Nektar::PreconCfs::m_CalcPreconMatFlag, Nektar::PreconCfs::m_Comm, Nektar::PreconCfs::m_DtLambdaPreconMat, Nektar::PreconCfs::m_operator, m_PreconItsStep, m_PreconMatVarsSingle, Nektar::PreconCfs::m_PreconTimesCounter, m_sBlkDiagMat, m_TraceIPSymJacArraySingle, m_TraceJacArraySingle, m_TraceJacDerivArraySingle, m_TraceJacDerivSignSingle, m_TraceJacDerivSingle, m_TraceJacSingle, and Nektar::PreconCfs::m_verbose.

◆ v_DoPreconCfs()

void Nektar::PreconCfsBRJ::v_DoPreconCfs	(	const Array< OneD, MultiRegions::ExpListSharedPtr > &	pFields,
		const Array< OneD, NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput,
		const bool &	flag
	)

overrideprotectedvirtual

Implements Nektar::PreconCfs.

Definition at line 85 of file PreconCfsBRJ.cpp.

{
    size_t nBRJIterTot = m_PreconItsStep;
    if (nBRJIterTot == 0)
    {
        DoNullPrecon(inarray, outarray, flag);
    }
    else
    {
        const NekDouble BRJParam   = m_BRJRelaxParam;
        const NekDouble OmBRJParam = 1.0 - BRJParam;
 
        size_t nvariables = pFields.size();
        size_t npoints    = pFields[0]->GetNcoeffs();
        size_t ntotpnt    = inarray.size();
 
        ASSERTL0(nvariables * npoints == ntotpnt,
                 "nvariables*npoints!=ntotpnt in PreconCoeff");
 
        Array<OneD, NekDouble> rhs(ntotpnt);
 
        Array<OneD, NekDouble> outN(ntotpnt);
        Array<OneD, NekDouble> outTmp(ntotpnt);
        Array<OneD, Array<OneD, NekDouble>> rhs2d(nvariables);
        Array<OneD, Array<OneD, NekDouble>> out_2d(nvariables);
        Array<OneD, Array<OneD, NekDouble>> outTmp_2d(nvariables);
        for (size_t m = 0; m < nvariables; m++)
        {
            size_t moffset = m * npoints;
            rhs2d[m]       = rhs + moffset;
            out_2d[m]      = outarray + moffset;
            outTmp_2d[m]   = outTmp + moffset;
            pFields[m]->MultiplyByMassMatrix(inarray + moffset, rhs2d[m]);
        }
 
        size_t nphysic   = pFields[0]->GetNpoints();
        size_t nTracePts = pFields[0]->GetTrace()->GetNpoints();
        TensorOfArray3D<NekDouble> qfield(m_spacedim);
        for (int i = 0; i < m_spacedim; i++)
        {
            qfield[i] = Array<OneD, Array<OneD, NekDouble>>(nvariables);
            for (size_t j = 0; j < nvariables; j++)
            {
                qfield[i][j] = Array<OneD, NekDouble>(nphysic);
            }
        }
        size_t ntmpTrace = 4 + 2 * m_spacedim;
        TensorOfArray3D<NekDouble> tmpTrace(ntmpTrace);
        for (size_t i = 0; i < ntmpTrace; i++)
        {
            tmpTrace[i] = Array<OneD, Array<OneD, NekDouble>>(nvariables);
            for (size_t j = 0; j < nvariables; j++)
            {
                tmpTrace[i][j] = Array<OneD, NekDouble>(nTracePts);
            }
        }
        Array<OneD, Array<OneD, NekDouble>> FwdFluxDeriv(nvariables);
        Array<OneD, Array<OneD, NekDouble>> BwdFluxDeriv(nvariables);
        for (size_t j = 0; j < nvariables; j++)
        {
            FwdFluxDeriv[j] = Array<OneD, NekDouble>(nTracePts);
            BwdFluxDeriv[j] = Array<OneD, NekDouble>(nTracePts);
        }
 
        const size_t nwspTraceDataType = nvariables + 1;
        Array<OneD, Array<OneD, NekSingle>> wspTraceDataType(nwspTraceDataType);
        for (size_t m = 0; m < nwspTraceDataType; m++)
        {
            wspTraceDataType[m] = Array<OneD, NekSingle>(nTracePts);
        }
 
        LibUtilities::Timer timer;
        timer.Start();
        PreconBlkDiag(pFields, rhs, outarray);
        timer.Stop();
        timer.AccumulateRegion("PreconCfsBRJ::PreconBlkDiag", 2);
 
        for (size_t nrelax = 0; nrelax < nBRJIterTot - 1; nrelax++)
        {
            Vmath::Smul(ntotpnt, OmBRJParam, outarray, 1, outN, 1);
 
            timer.Start();
            MinusOffDiag2Rhs(pFields, nvariables, npoints, rhs2d, out_2d,
                             tmpTrace, wspTraceDataType, m_TraceJacArraySingle);
            timer.Stop();
            timer.AccumulateRegion("PreconCfsBRJ::MinusOffDiag2Rhs", 2);
 
            timer.Start();
            PreconBlkDiag(pFields, outarray, outTmp);
            timer.Stop();
            timer.AccumulateRegion("PreconCfsBRJ::PreconBlkDiag", 2);
 
            Vmath::Svtvp(ntotpnt, BRJParam, outTmp, 1, outN, 1, outarray, 1);
        }
    }
}

References Nektar::LibUtilities::Timer::AccumulateRegion(), ASSERTL0, DoNullPrecon(), m_BRJRelaxParam, m_PreconItsStep, Nektar::PreconCfs::m_spacedim, m_TraceJacArraySingle, MinusOffDiag2Rhs(), PreconBlkDiag(), Vmath::Smul(), Nektar::LibUtilities::Timer::Start(), Nektar::LibUtilities::Timer::Stop(), and Vmath::Svtvp().

◆ v_InitObject()

void Nektar::PreconCfsBRJ::v_InitObject ( )

overrideprotectedvirtual

Implements Nektar::PreconCfs.

Definition at line 78 of file PreconCfsBRJ.cpp.

79{

80}

◆ v_UpdatePreconMatCheck()

bool Nektar::PreconCfsBRJ::v_UpdatePreconMatCheck	(	const Array< OneD, const NekDouble > &	res,
		const NekDouble	dtLambda
	)

overrideprotectedvirtual

Implements Nektar::PreconCfs.

Definition at line 275 of file PreconCfsBRJ.cpp.

{
    bool flag = (m_CalcPreconMatFlag || m_DtLambdaPreconMat != dtLambda ||
                 m_PreconMatFreezNumb < m_PreconTimesCounter);
    m_CalcPreconMatFlag = flag;
    return flag;
}

References Nektar::PreconCfs::m_CalcPreconMatFlag, Nektar::PreconCfs::m_DtLambdaPreconMat, Nektar::PreconCfs::m_PreconMatFreezNumb, and Nektar::PreconCfs::m_PreconTimesCounter.

Friends And Related Function Documentation

◆ MemoryManager< PreconCfsBRJ >

friend class MemoryManager< PreconCfsBRJ >

friend

Definition at line 1 of file PreconCfsBRJ.h.

Member Data Documentation

◆ className

std::string Nektar::PreconCfsBRJ::className

static

Initial value:

=
    GetPreconCfsFactory().RegisterCreatorFunction(
        "PreconCfsBRJ", PreconCfsBRJ::create,
        "Block Relaxed Jacobi Preconditioner for CFS.")

Name of the class.

Definition at line 66 of file PreconCfsBRJ.h.

◆ m_BRJRelaxParam

int Nektar::PreconCfsBRJ::m_BRJRelaxParam

protected

Definition at line 75 of file PreconCfsBRJ.h.

Referenced by PreconCfsBRJ(), and v_DoPreconCfs().

◆ m_inputIdx

std::vector<int> Nektar::PreconCfsBRJ::m_inputIdx

protected

Definition at line 83 of file PreconCfsBRJ.h.

Referenced by AllocateSIMDPreconBlkMatDiag(), and PreconBlkDiag().

◆ m_max_nblocks

unsigned int Nektar::PreconCfsBRJ::m_max_nblocks

protected

Definition at line 79 of file PreconCfsBRJ.h.

Referenced by AllocateSIMDPreconBlkMatDiag(), and PreconBlkDiag().

◆ m_max_nElmtDof

unsigned int Nektar::PreconCfsBRJ::m_max_nElmtDof

protected

Definition at line 80 of file PreconCfsBRJ.h.

Referenced by AllocateSIMDPreconBlkMatDiag(), and PreconBlkDiag().

◆ m_PreconItsStep

int Nektar::PreconCfsBRJ::m_PreconItsStep

protected

Definition at line 74 of file PreconCfsBRJ.h.

Referenced by PreconCfsBRJ(), v_BuildPreconCfs(), and v_DoPreconCfs().

◆ m_PreconMatVarsSingle

Array<OneD, Array<OneD, SNekBlkMatSharedPtr> > Nektar::PreconCfsBRJ::m_PreconMatVarsSingle

protected

Definition at line 77 of file PreconCfsBRJ.h.

Referenced by PreconCfsBRJ(), and v_BuildPreconCfs().

◆ m_sBlkDiagMat

std::vector<simd<NekSingle>, tinysimd::allocator<simd<NekSingle> > > Nektar::PreconCfsBRJ::m_sBlkDiagMat

protected

Definition at line 82 of file PreconCfsBRJ.h.

Referenced by AllocateSIMDPreconBlkMatDiag(), PreconBlkDiag(), and v_BuildPreconCfs().

◆ m_TraceIPSymJacArraySingle

TensorOfArray5D<NekSingle> Nektar::PreconCfsBRJ::m_TraceIPSymJacArraySingle

protected

Definition at line 90 of file PreconCfsBRJ.h.

Referenced by v_BuildPreconCfs().

◆ m_TraceJacArraySingle

TensorOfArray4D<NekSingle> Nektar::PreconCfsBRJ::m_TraceJacArraySingle

protected

Definition at line 86 of file PreconCfsBRJ.h.

Referenced by v_BuildPreconCfs(), and v_DoPreconCfs().

◆ m_TraceJacDerivArraySingle

TensorOfArray4D<NekSingle> Nektar::PreconCfsBRJ::m_TraceJacDerivArraySingle

protected

Definition at line 88 of file PreconCfsBRJ.h.

Referenced by v_BuildPreconCfs().

◆ m_TraceJacDerivSignSingle

Array<OneD, Array<OneD, NekSingle> > Nektar::PreconCfsBRJ::m_TraceJacDerivSignSingle

protected

Definition at line 89 of file PreconCfsBRJ.h.

Referenced by v_BuildPreconCfs().

◆ m_TraceJacDerivSingle

Array<OneD, SNekBlkMatSharedPtr> Nektar::PreconCfsBRJ::m_TraceJacDerivSingle

protected

Definition at line 87 of file PreconCfsBRJ.h.

Referenced by v_BuildPreconCfs().

◆ m_TraceJacSingle

Array<OneD, SNekBlkMatSharedPtr> Nektar::PreconCfsBRJ::m_TraceJacSingle

protected

Definition at line 85 of file PreconCfsBRJ.h.

Referenced by v_BuildPreconCfs().

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Protected Attributes

Private Member Functions

Friends

Detailed Description

Constructor & Destructor Documentation

◆ PreconCfsBRJ()

◆ ~PreconCfsBRJ()

Member Function Documentation

◆ AllocateNekBlkMatDig()

◆ AllocatePreconBlkDiagCoeff()

◆ AllocateSIMDPreconBlkMatDiag()

◆ create()

◆ DoNullPrecon()

◆ MinusOffDiag2Rhs()

◆ PreconBlkDiag()

◆ v_BuildPreconCfs()

◆ v_DoPreconCfs()

◆ v_InitObject()

◆ v_UpdatePreconMatCheck()

Friends And Related Function Documentation

◆ MemoryManager< PreconCfsBRJ >

Member Data Documentation

◆ className

◆ m_BRJRelaxParam

◆ m_inputIdx

◆ m_max_nblocks

◆ m_max_nElmtDof

◆ m_PreconItsStep

◆ m_PreconMatVarsSingle

◆ m_sBlkDiagMat

◆ m_TraceIPSymJacArraySingle

◆ m_TraceJacArraySingle

◆ m_TraceJacDerivArraySingle

◆ m_TraceJacDerivSignSingle

◆ m_TraceJacDerivSingle

◆ m_TraceJacSingle