A global linear system. More...

#include <GlobalLinSysIterativeFull.h>

Inheritance diagram for Nektar::MultiRegions::GlobalLinSysIterativeFull:

Public Member Functions
	GlobalLinSysIterativeFull (const GlobalLinSysKey &pLinSysKey, const std::weak_ptr< ExpList > &pExpList, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
	Constructor for full direct matrix solve. More...

	~GlobalLinSysIterativeFull () override=default

Public Member Functions inherited from Nektar::MultiRegions::GlobalLinSysIterative
	GlobalLinSysIterative (const GlobalLinSysKey &pKey, const std::weak_ptr< ExpList > &pExpList, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
	Constructor for full direct matrix solve. More...

	~GlobalLinSysIterative () override

void	DoMatrixMultiply (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)

Public Member Functions inherited from Nektar::MultiRegions::GlobalLinSys
	GlobalLinSys (const GlobalLinSysKey &pKey, const std::weak_ptr< ExpList > &pExpList, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
	Constructor for full direct matrix solve. More...

virtual	~GlobalLinSys ()

const GlobalLinSysKey &	GetKey (void) const
	Returns the key associated with the system. More...

const std::weak_ptr< ExpList > &	GetLocMat (void) const

void	InitObject ()

void	Initialise (const std::shared_ptr< AssemblyMap > &pLocToGloMap)

void	Solve (const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out, const AssemblyMapSharedPtr &locToGloMap, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
	Solve the linear system for given input and output vectors using a specified local to global map. More...

std::shared_ptr< GlobalLinSys >	GetSharedThisPtr ()
	Returns a shared pointer to the current object. More...

int	GetNumBlocks ()

DNekScalMatSharedPtr	GetBlock (unsigned int n)

void	DropBlock (unsigned int n)

DNekScalBlkMatSharedPtr	GetStaticCondBlock (unsigned int n)

void	DropStaticCondBlock (unsigned int n)

void	SolveLinearSystem (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir=0)
	Solve the linear system for given input and output vectors. More...

Static Public Member Functions
static GlobalLinSysSharedPtr	create (const GlobalLinSysKey &pLinSysKey, const std::weak_ptr< ExpList > &pExpList, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
	Creates an instance of this class. More...

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

Protected Member Functions
void	v_Solve (const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out, const AssemblyMapSharedPtr &locToGloMap, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray) override
	Solve the linear system for given input and output vectors using a specified local to global map. More...

void	v_DoMatrixMultiply (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput) override

void	v_UniqueMap () override

Protected Member Functions inherited from Nektar::MultiRegions::GlobalLinSysIterative
void	DoProjection (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const int pNumDir, const NekDouble tol, const bool isAconjugate)
	projection technique More...

void	Set_Rhs_Magnitude (const NekVector< NekDouble > &pIn)

virtual void	v_UniqueMap ()=0

virtual void	v_DoMatrixMultiply (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)=0

Protected Member Functions inherited from Nektar::MultiRegions::GlobalLinSys
virtual void	v_Solve (const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out, const AssemblyMapSharedPtr &locToGloMap, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)=0
	Solve a linear system based on mapping. More...

virtual void	v_SolveLinearSystem (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)=0
	Solve a basic matrix system. More...

virtual void	v_InitObject ()

virtual void	v_Initialise (const std::shared_ptr< AssemblyMap > &pLocToGloMap)

virtual int	v_GetNumBlocks ()
	Get the number of blocks in this system. More...

virtual DNekScalMatSharedPtr	v_GetBlock (unsigned int n)
	Retrieves the block matrix from n-th expansion using the matrix key provided by the m_linSysKey. More...

virtual void	v_DropBlock (unsigned int n)
	Releases the local block matrix from NekManager of n-th expansion using the matrix key provided by the m_linSysKey. More...

virtual DNekScalBlkMatSharedPtr	v_GetStaticCondBlock (unsigned int n)
	Retrieves a the static condensation block matrices from n-th expansion using the matrix key provided by the m_linSysKey. More...

virtual void	v_DropStaticCondBlock (unsigned int n)
	Releases the static condensation block matrices from NekManager of n-th expansion using the matrix key provided by the m_linSysKey. More...

PreconditionerSharedPtr	CreatePrecon (AssemblyMapSharedPtr asmMap)
	Create a preconditioner object from the parameters defined in the supplied assembly map. More...

Private Member Functions
void	v_SolveLinearSystem (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir) override
	Solve the matrix system. More...

Private Attributes
std::weak_ptr< AssemblyMap >	m_locToGloMap

Additional Inherited Members
Protected Attributes inherited from Nektar::MultiRegions::GlobalLinSysIterative
Array< OneD, int >	m_map
	Global to universal unique map. More...

NekDouble	m_tolerance
	Tolerance of iterative solver. More...

NekDouble	m_rhs_magnitude
	dot product of rhs to normalise stopping criterion More...

NekDouble	m_rhs_mag_sm
	cnt to how many times rhs_magnitude is called More...

PreconditionerSharedPtr	m_precon

std::string	m_precontype

int	m_totalIterations

bool	m_useProjection
	Whether to apply projection technique. More...

bool	m_root
	Root if parallel. More...

std::string	m_linSysIterSolver
	Iterative solver: Conjugate Gradient, GMRES. More...

std::vector< Array< OneD, NekDouble > >	m_prevLinSol
	Storage for solutions to previous linear problems. More...

std::vector< Array< OneD, NekDouble > >	m_prevBasis

DNekMatSharedPtr	m_coeffMatrix

Array< OneD, NekDouble >	m_coeffMatrixFactor

Array< OneD, int >	m_ipivot

int	m_numSuccessiveRHS

bool	m_isAconjugate

std::string	m_matrixType

bool	m_isNonSymmetricLinSys

int	m_numPrevSols

bool	m_isAbsoluteTolerance

LibUtilities::NekSysOperators	m_NekSysOp

LibUtilities::NekLinSysIterSharedPtr	m_linsol

Protected Attributes inherited from Nektar::MultiRegions::GlobalLinSys
const GlobalLinSysKey	m_linSysKey
	Key associated with this linear system. More...

const std::weak_ptr< ExpList >	m_expList
	Local Matrix System. More...

const std::map< int, RobinBCInfoSharedPtr >	m_robinBCInfo
	Robin boundary info. More...

bool	m_verbose

Static Protected Attributes inherited from Nektar::MultiRegions::GlobalLinSysIterative
static std::string	IteratSolverlookupIds []

static std::string	IteratSolverdef

Detailed Description

A global linear system.

Definition at line 45 of file GlobalLinSysIterativeFull.h.

Constructor & Destructor Documentation

◆ GlobalLinSysIterativeFull()

Nektar::MultiRegions::GlobalLinSysIterativeFull::GlobalLinSysIterativeFull	(	const GlobalLinSysKey &	pKey,
		const std::weak_ptr< ExpList > &	pExp,
		const std::shared_ptr< AssemblyMap > &	pLocToGloMap
	)

Constructor for full direct matrix solve.

Parameters

pKey	Key specifying matrix to solve.
pExp	Shared pointer to expansion list for applying matrix evaluations.
pLocToGloMap	Local to global mapping.

Definition at line 66 of file GlobalLinSysIterativeFull.cpp.

    : GlobalLinSys(pKey, pExp, pLocToGloMap),
      GlobalLinSysIterative(pKey, pExp, pLocToGloMap)
{
    ASSERTL1(m_linSysKey.GetGlobalSysSolnType() == eIterativeFull,
             "This routine should only be used when using an Iterative "
             "conjugate gradient matrix solve.");
}

References ASSERTL1, Nektar::MultiRegions::eIterativeFull, Nektar::MultiRegions::GlobalLinSysKey::GetGlobalSysSolnType(), and Nektar::MultiRegions::GlobalLinSys::m_linSysKey.

◆ ~GlobalLinSysIterativeFull()

Nektar::MultiRegions::GlobalLinSysIterativeFull::~GlobalLinSysIterativeFull ( )

overridedefault

Member Function Documentation

◆ create()

static GlobalLinSysSharedPtr Nektar::MultiRegions::GlobalLinSysIterativeFull::create	(	const GlobalLinSysKey &	pLinSysKey,
		const std::weak_ptr< ExpList > &	pExpList,
		const std::shared_ptr< AssemblyMap > &	pLocToGloMap
	)

inlinestatic

Creates an instance of this class.

Definition at line 49 of file GlobalLinSysIterativeFull.h.

    {
        return MemoryManager<GlobalLinSysIterativeFull>::AllocateSharedPtr(
            pLinSysKey, pExpList, pLocToGloMap);
    }

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

◆ v_DoMatrixMultiply()

void Nektar::MultiRegions::GlobalLinSysIterativeFull::v_DoMatrixMultiply	(	const Array< OneD, NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput
	)

overrideprotectedvirtual

Implements Nektar::MultiRegions::GlobalLinSysIterative.

Definition at line 179 of file GlobalLinSysIterativeFull.cpp.

{
    bool isLocal = m_linsol->IsLocal();
 
    std::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
 
    AssemblyMapSharedPtr asmMap = m_locToGloMap.lock();
 
    int ncoeffs = expList->GetNcoeffs();
    Array<OneD, NekDouble> InputLoc, OutputLoc;
 
    if (isLocal)
    {
        InputLoc  = pInput;
        OutputLoc = pOutput;
 
        // Perform matrix-vector operation A*d_i
        expList->GeneralMatrixOp(m_linSysKey, InputLoc, OutputLoc);
    }
    else
    {
        InputLoc  = Array<OneD, NekDouble>(ncoeffs);
        OutputLoc = Array<OneD, NekDouble>(ncoeffs);
 
        asmMap->GlobalToLocal(pInput, InputLoc);
        // Perform matrix-vector operation A*d_i
        expList->GeneralMatrixOp(m_linSysKey, InputLoc, OutputLoc);
    }
 
    // Apply robin boundary conditions to the solution.
    for (auto &r : m_robinBCInfo) // add robin mass matrix
    {
        RobinBCInfoSharedPtr rBC;
        Array<OneD, NekDouble> tmp;
 
        int n = r.first;
 
        int offset                            = expList->GetCoeff_Offset(n);
        LocalRegions::ExpansionSharedPtr vExp = expList->GetExp(n);
 
        // add local matrix contribution
        for (rBC = r.second; rBC; rBC = rBC->next)
        {
            vExp->AddRobinTraceContribution(
                rBC->m_robinID, rBC->m_robinPrimitiveCoeffs, InputLoc + offset,
                tmp = OutputLoc + offset);
        }
    }
 
    // put back in global coeffs
    if (isLocal == false)
    {
        asmMap->Assemble(OutputLoc, pOutput);
    }
}

References Nektar::MultiRegions::GlobalLinSys::m_expList, Nektar::MultiRegions::GlobalLinSysIterative::m_linsol, Nektar::MultiRegions::GlobalLinSys::m_linSysKey, m_locToGloMap, and Nektar::MultiRegions::GlobalLinSys::m_robinBCInfo.

◆ v_Solve()

void Nektar::MultiRegions::GlobalLinSysIterativeFull::v_Solve	(	const Array< OneD, const NekDouble > &	pLocInput,
		Array< OneD, NekDouble > &	pLocOutput,
		const AssemblyMapSharedPtr &	pLocToGloMap,
		const Array< OneD, const NekDouble > &	pDirForcing = `NullNekDouble1DArray`
	)

overrideprotectedvirtual

Solve the linear system for given input and output vectors using a specified local to global map.

Solve a global linear system with Dirichlet forcing using a conjugate gradient method. This routine performs handling of the Dirichlet forcing terms and wraps the underlying iterative solver used for the remaining degrees of freedom.

Consider solving for \(x\), the matrix system \(Ax=b\), where \(b\) is known. To enforce the Dirichlet terms we instead solve

\[A(x-x_0) = b - Ax_0 \]

where \(x_0\) is the Dirichlet forcing.

Parameters

pInput	RHS of linear system, \(b\).
pOutput	On input, values of dirichlet degrees of freedom with initial guess on other values. On output, the solution \(x\).
pLocToGloMap	Local to global mapping.
pDirForcing	Precalculated Dirichlet forcing.

Implements Nektar::MultiRegions::GlobalLinSys.

Definition at line 95 of file GlobalLinSysIterativeFull.cpp.

{
    m_locToGloMap = pLocToGloMap;
 
    bool dirForcCalculated = (bool)pDirForcing.size();
    int nDirDofs           = pLocToGloMap->GetNumGlobalDirBndCoeffs();
    int nGlobDofs          = pLocToGloMap->GetNumGlobalCoeffs();
    int nLocDofs           = pLocToGloMap->GetNumLocalCoeffs();
 
    int nDirTotal                                  = nDirDofs;
    std::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
    expList->GetComm()->GetRowComm()->AllReduce(nDirTotal,
                                                LibUtilities::ReduceSum);
 
    if (nDirTotal)
    {
        Array<OneD, NekDouble> rhs(nLocDofs);
 
        // Calculate the Dirichlet forcing
        if (dirForcCalculated)
        {
            // Assume pDirForcing is in local space
            ASSERTL0(
                pDirForcing.size() >= nLocDofs,
                "DirForcing is not of sufficient size. Is it in local space?");
            Vmath::Vsub(nLocDofs, pLocInput, 1, pDirForcing, 1, rhs, 1);
        }
        else
        {
            // Calculate initial condition and Dirichlet forcing and subtract it
            // from the rhs
            expList->GeneralMatrixOp(m_linSysKey, pLocOutput, rhs);
 
            // Iterate over all the elements computing Robin BCs where
            // necessary
            for (auto &r : m_robinBCInfo) // add robin mass matrix
            {
                RobinBCInfoSharedPtr rBC;
                Array<OneD, NekDouble> rhsloc;
 
                int n      = r.first;
                int offset = expList->GetCoeff_Offset(n);
 
                LocalRegions::ExpansionSharedPtr vExp = expList->GetExp(n);
                // Add local matrix contribution
                for (rBC = r.second; rBC; rBC = rBC->next)
                {
                    vExp->AddRobinTraceContribution(
                        rBC->m_robinID, rBC->m_robinPrimitiveCoeffs,
                        pLocOutput + offset, rhsloc = rhs + offset);
                }
            }
            Vmath::Vsub(nLocDofs, pLocInput, 1, rhs, 1, rhs, 1);
        }
 
        if (std::dynamic_pointer_cast<AssemblyMapCG>(pLocToGloMap))
        {
            Array<OneD, NekDouble> diff(nLocDofs);
 
            // Solve for perturbation from initial guess in pOutput
            SolveLinearSystem(nGlobDofs, rhs, diff, pLocToGloMap, nDirDofs);
 
            // Add back initial and boundary condition
            Vmath::Vadd(nLocDofs, diff, 1, pLocOutput, 1, pLocOutput, 1);
        }
        else
        {
            ASSERTL0(false, "Need DG solve if using Dir BCs");
        }
    }
    else
    {
        SolveLinearSystem(nGlobDofs, pLocInput, pLocOutput, pLocToGloMap,
                          nDirDofs);
    }
}

References ASSERTL0, Nektar::MultiRegions::GlobalLinSys::m_expList, Nektar::MultiRegions::GlobalLinSys::m_linSysKey, m_locToGloMap, Nektar::MultiRegions::GlobalLinSys::m_robinBCInfo, Nektar::LibUtilities::ReduceSum, Nektar::MultiRegions::GlobalLinSys::SolveLinearSystem(), Vmath::Vadd(), and Vmath::Vsub().

◆ v_SolveLinearSystem()

void Nektar::MultiRegions::GlobalLinSysIterativeFull::v_SolveLinearSystem	(	const int	pNumRows,
		const Array< OneD, const NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput,
		const AssemblyMapSharedPtr &	locToGloMap,
		const int	pNumDir
	)

overrideprivatevirtual

Solve the matrix system.

Implements Nektar::MultiRegions::GlobalLinSys.

Definition at line 247 of file GlobalLinSysIterativeFull.cpp.

{
 
    if (!m_linsol)
    {
        LibUtilities::CommSharedPtr vRowComm =
            m_expList.lock()->GetComm()->GetRowComm();
        LibUtilities::SessionReaderSharedPtr pSession =
            m_expList.lock()->GetSession();
 
        // Check such a module exists for this equation.
        ASSERTL0(LibUtilities::GetNekLinSysIterFactory().ModuleExists(
                     m_linSysIterSolver),
                 "NekLinSysIter '" + m_linSysIterSolver +
                     "' is not defined.\n");
 
        // Create the key to hold solver settings
        auto sysKey     = LibUtilities::NekSysKey();
        string variable = plocToGloMap->GetVariable();
 
        // Either get the solnInfo from <GlobalSysSolInfo> or from
        // <Parameters>
        if (pSession->DefinesGlobalSysSolnInfo(variable,
                                               "NekLinSysMaxIterations"))
        {
            sysKey.m_NekLinSysMaxIterations = boost::lexical_cast<int>(
                pSession
                    ->GetGlobalSysSolnInfo(variable, "NekLinSysMaxIterations")
                    .c_str());
        }
        else
        {
            pSession->LoadParameter("NekLinSysMaxIterations",
                                    sysKey.m_NekLinSysMaxIterations, 5000);
        }
 
        if (pSession->DefinesGlobalSysSolnInfo(variable, "LinSysMaxStorage"))
        {
            sysKey.m_LinSysMaxStorage = boost::lexical_cast<int>(
                pSession->GetGlobalSysSolnInfo(variable, "LinSysMaxStorage")
                    .c_str());
        }
        else
        {
            pSession->LoadParameter("LinSysMaxStorage",
                                    sysKey.m_LinSysMaxStorage, 100);
        }
 
        if (pSession->DefinesGlobalSysSolnInfo(variable, "GMRESMaxHessMatBand"))
        {
            sysKey.m_KrylovMaxHessMatBand = boost::lexical_cast<int>(
                pSession->GetGlobalSysSolnInfo(variable, "GMRESMaxHessMatBand")
                    .c_str());
        }
        else
        {
            pSession->LoadParameter("GMRESMaxHessMatBand",
                                    sysKey.m_KrylovMaxHessMatBand,
                                    sysKey.m_LinSysMaxStorage + 1);
        }
 
        // The following settings have no correponding tests and are rarely
        // used.
        pSession->MatchSolverInfo("GMRESLeftPrecon", "True",
                                  sysKey.m_NekLinSysLeftPrecon, false);
        pSession->MatchSolverInfo("GMRESRightPrecon", "True",
                                  sysKey.m_NekLinSysRightPrecon, true);
 
        m_linsol = LibUtilities::GetNekLinSysIterFactory().CreateInstance(
            m_linSysIterSolver, pSession, vRowComm, nGlobal - nDir, sysKey);
 
        m_linsol->SetSysOperators(m_NekSysOp);
        v_UniqueMap();
        m_linsol->SetUniversalUniqueMap(m_map);
    }
 
    if (!m_precon)
    {
        m_precon = CreatePrecon(plocToGloMap);
        m_precon->BuildPreconditioner();
    }
 
    m_linsol->setRhsMagnitude(m_isAbsoluteTolerance ? 1.0 : m_rhs_magnitude);
 
    if (m_useProjection)
    {
        Array<OneD, NekDouble> gloIn(nGlobal);
        Array<OneD, NekDouble> gloOut(nGlobal, 0.0);
        plocToGloMap->Assemble(pInput, gloIn);
        DoProjection(nGlobal, gloIn, gloOut, nDir, m_tolerance, m_isAconjugate);
        plocToGloMap->GlobalToLocal(gloOut, pOutput);
    }
    else
    {
        if (m_linsol->IsLocal())
        {
            int nLocDofs = plocToGloMap->GetNumLocalCoeffs();
            Vmath::Zero(nLocDofs, pOutput, 1);
            m_linsol->SolveSystem(nLocDofs, pInput, pOutput, nDir, m_tolerance);
        }
        else
        {
            Array<OneD, NekDouble> gloIn(nGlobal);
            Array<OneD, NekDouble> gloOut(nGlobal, 0.0);
            plocToGloMap->Assemble(pInput, gloIn);
            m_linsol->SolveSystem(nGlobal, gloIn, gloOut, nDir, m_tolerance);
            plocToGloMap->GlobalToLocal(gloOut, pOutput);
        }
    }
}

◆ v_UniqueMap()

void Nektar::MultiRegions::GlobalLinSysIterativeFull::v_UniqueMap ( )

overrideprotectedvirtual

Implements Nektar::MultiRegions::GlobalLinSysIterative.

Definition at line 239 of file GlobalLinSysIterativeFull.cpp.

{
    m_map = m_locToGloMap.lock()->GetGlobalToUniversalMapUnique();
}

References m_locToGloMap, and Nektar::MultiRegions::GlobalLinSysIterative::m_map.

Referenced by v_SolveLinearSystem().

Member Data Documentation

◆ className

string Nektar::MultiRegions::GlobalLinSysIterativeFull::className

static

Initial value:

=
    GetGlobalLinSysFactory().RegisterCreatorFunction(
        "IterativeFull", GlobalLinSysIterativeFull::create,
        "Iterative solver for full matrix system.")

Name of class.

Registers the class with the Factory.

Definition at line 59 of file GlobalLinSysIterativeFull.h.

◆ m_locToGloMap

std::weak_ptr<AssemblyMap> Nektar::MultiRegions::GlobalLinSysIterativeFull::m_locToGloMap

private

Definition at line 85 of file GlobalLinSysIterativeFull.h.

Referenced by v_DoMatrixMultiply(), v_Solve(), and v_UniqueMap().

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ GlobalLinSysIterativeFull()

◆ ~GlobalLinSysIterativeFull()

Member Function Documentation

◆ create()

◆ v_DoMatrixMultiply()

◆ v_Solve()

◆ v_SolveLinearSystem()

◆ v_UniqueMap()

Member Data Documentation

◆ className

◆ m_locToGloMap