Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Attributes | Private Member Functions | List of all members
Nektar::MultiRegions::PreconditionerBlock Class Reference

#include <PreconditionerBlock.h>

Inheritance diagram for Nektar::MultiRegions::PreconditionerBlock:
Inheritance graph
[legend]
Collaboration diagram for Nektar::MultiRegions::PreconditionerBlock:
Collaboration graph
[legend]

Public Member Functions

 PreconditionerBlock (const boost::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
virtual ~PreconditionerBlock ()
- Public Member Functions inherited from Nektar::MultiRegions::Preconditioner
 Preconditioner (const boost::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
virtual ~Preconditioner ()
void DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
void DoPreconditionerWithNonVertOutput (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const Array< OneD, NekDouble > &pNonVertOutput, Array< OneD, NekDouble > &pVertForce=NullNekDouble1DArray)
void DoTransformToLowEnergy (Array< OneD, NekDouble > &pInOut, int offset)
void DoTransformToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
void DoTransformFromLowEnergy (Array< OneD, NekDouble > &pInOut)
void DoMultiplybyInverseTransformationMatrix (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
void DoMultiplybyInverseTransposedTransformationMatrix (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
void BuildPreconditioner ()
void InitObject ()
Array< OneD, NekDouble > AssembleStaticCondGlobalDiagonals ()
 Performs global assembly of diagonal entries to global Schur complement matrix.
const DNekScalBlkMatSharedPtr & GetBlockTransformedSchurCompl () const
const DNekScalBlkMatSharedPtr & GetBlockCMatrix () const
const DNekScalBlkMatSharedPtr & GetBlockInvDMatrix () const
const DNekScalBlkMatSharedPtr & GetBlockSchurCompl () const
const DNekScalBlkMatSharedPtr & GetBlockTransformationMatrix () const
const DNekScalBlkMatSharedPtr & GetBlockTransposedTransformationMatrix () const
DNekScalMatSharedPtr TransformedSchurCompl (int offset, const boost::shared_ptr< DNekScalMat > &loc_mat)

Static Public Member Functions

static PreconditionerSharedPtr create (const boost::shared_ptr< GlobalLinSys > &plinsys, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
 Creates an instance of this class.

Static Public Attributes

static std::string className
 Name of class.

Protected Attributes

const boost::weak_ptr
< GlobalLinSys > 		m_linsys
PreconditionerType m_preconType
DNekBlkMatSharedPtr m_blkMat
boost::shared_ptr< AssemblyMap > m_locToGloMap
- Protected Attributes inherited from Nektar::MultiRegions::Preconditioner
const boost::weak_ptr
< GlobalLinSys > 		m_linsys
PreconditionerType m_preconType
DNekMatSharedPtr m_preconditioner
boost::shared_ptr< AssemblyMap > m_locToGloMap
LibUtilities::CommSharedPtr m_comm

Private Member Functions

void BlockPreconditioner2D (void)
 Construct a block preconditioner from $\mathbf{S}_{1}$.
void BlockPreconditioner3D (void)
void BlockPreconditionerHDG (void)
virtual void v_InitObject ()
virtual void v_DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
virtual void v_BuildPreconditioner ()

Additional Inherited Members

- Protected Member Functions inherited from Nektar::MultiRegions::Preconditioner
virtual DNekScalMatSharedPtr v_TransformedSchurCompl (int offset, const boost::shared_ptr< DNekScalMat > &loc_mat)
 Get block elemental transposed transformation matrix $\mathbf{R}^{T}$.

Detailed Description

Definition at line 53 of file PreconditionerBlock.h.

Constructor & Destructor Documentation

Nektar::MultiRegions::PreconditionerBlock::PreconditionerBlock ( const boost::shared_ptr< GlobalLinSys > &  plinsys,
const AssemblyMapSharedPtr &  pLocToGloMap 
)

Definition at line 65 of file PreconditionerBlock.cpp.

: Preconditioner(plinsys, pLocToGloMap),
m_linsys(plinsys),
m_preconType(pLocToGloMap->GetPreconType()),
m_locToGloMap(pLocToGloMap)
{
}
virtual Nektar::MultiRegions::PreconditionerBlock::~PreconditionerBlock ( )
inlinevirtual

Definition at line 75 of file PreconditionerBlock.h.

{}

Member Function Documentation

void Nektar::MultiRegions::PreconditionerBlock::BlockPreconditioner2D ( void  )
private

Construct a block preconditioner from $\mathbf{S}_{1}$.

\[\mathbf{M}^{-1}=\left[\begin{array}{ccc} Diag[(\mathbf{S_{1}})_{vv}] & & \\ & (\mathbf{S}_{1})_{eb} & \\ & & (\mathbf{S}_{1})_{fb} \end{array}\right] \]

where $\mathbf{S}_{1}$ is the local schur complement matrix for each element.

Definition at line 123 of file PreconditionerBlock.cpp.

References Vmath::Assmb(), Nektar::eDIAGONAL, Nektar::SpatialDomains::eDirichlet, Nektar::eFULL, Gs::Gather(), Nektar::StdRegions::StdExpansion::GetEdgeInverseBoundaryMap(), Nektar::StdRegions::StdExpansion::GetEdgeNcoeffs(), Gs::gs_add, Gs::Init(), m_blkMat, Nektar::MultiRegions::Preconditioner::m_comm, m_linsys, m_locToGloMap, Nektar::LibUtilities::ReduceMax, and Vmath::Zero().

Referenced by v_BuildPreconditioner().

{
boost::shared_ptr<MultiRegions::ExpList>
expList=((m_linsys.lock())->GetLocMat()).lock();
LocalRegions::ExpansionSharedPtr locExpansion;
GlobalLinSysKey m_linSysKey=(m_linsys.lock())->GetKey();
StdRegions::VarCoeffMap vVarCoeffMap;
int i, j, k, nel;
int nVerts, nEdges;
int eid, n, cnt, nedgemodes;
NekDouble zero = 0.0;
int vMap1, vMap2, sign1, sign2;
int m, v, eMap1, eMap2;
int offset, globalrow, globalcol;
//matrix storage
MatrixStorage storage = eFULL;
MatrixStorage vertstorage = eDIAGONAL;
MatrixStorage blkmatStorage = eDIAGONAL;
//local element static condensed matrices
DNekScalBlkMatSharedPtr loc_mat;
DNekScalMatSharedPtr bnd_mat;
DNekMatSharedPtr VertBlk;
int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
int nNonDirVerts = m_locToGloMap->GetNumNonDirVertexModes();
//Vertex and edge preconditioner matrices
VertBlk = MemoryManager<DNekMat>::
AllocateSharedPtr(nNonDirVerts,nNonDirVerts,zero,vertstorage);
Array<OneD, NekDouble> vertArray(nNonDirVerts,0.0);
Array<OneD, long> VertBlockToUniversalMap(nNonDirVerts,-1);
int n_exp = expList->GetNumElmts();
int nNonDirEdgeIDs=m_locToGloMap->GetNumNonDirEdges();
//set the number of blocks in the matrix
Array<OneD,unsigned int> n_blks(1+nNonDirEdgeIDs);
n_blks[0]=nNonDirVerts;
map<int,int> edgeDirMap;
map<int,int> uniqueEdgeMap;
//this should be of size total number of local edges
Array<OneD, int> edgemodeoffset(nNonDirEdgeIDs,0);
Array<OneD, int> edgeglobaloffset(nNonDirEdgeIDs,0);
const Array<OneD, const ExpListSharedPtr>& bndCondExp = expList->GetBndCondExpansions();
StdRegions::StdExpansion1DSharedPtr bndCondFaceExp;
LocalRegions::SegExpSharedPtr bndSegExp;
const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>&
bndConditions = expList->GetBndConditions();
int meshVertId;
int meshEdgeId;
//Determine which boundary edges and faces have dirichlet values
for(i = 0; i < bndCondExp.num_elements(); i++)
{
cnt = 0;
for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
{
bndSegExp = bndCondExp[i]->GetExp(j)
->as<LocalRegions::SegExp>();
if (bndConditions[i]->GetBoundaryConditionType() ==
SpatialDomains::eDirichlet)
{
meshEdgeId = (bndSegExp->GetGeom1D())->GetEid();
edgeDirMap[meshEdgeId] = 1;
}
}
}
int dof=0;
int maxEdgeDof=0;
int nlocalNonDirEdges=0;
int edgematrixlocation=0;
int ntotaledgeentries=0;
// Loop over all the elements in the domain and compute max edge
// DOF. Reduce across all processes to get universal maximum.
for(cnt=n=0; n < n_exp; ++n)
{
nel = expList->GetOffset_Elmt_Id(n);
locExpansion = expList->GetExp(nel);
for (j = 0; j < locExpansion->GetNedges(); ++j)
{
dof = locExpansion->GetEdgeNcoeffs(j)-2;
maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
meshEdgeId = locExpansion->as<LocalRegions::Expansion2D>()
->GetGeom2D()->GetEid(j);
if(edgeDirMap.count(meshEdgeId)==0)
{
if(uniqueEdgeMap.count(meshEdgeId)==0)
{
uniqueEdgeMap[meshEdgeId]=edgematrixlocation;
edgeglobaloffset[edgematrixlocation]+=ntotaledgeentries;
edgemodeoffset[edgematrixlocation]=dof*dof;
ntotaledgeentries+=dof*dof;
n_blks[1+edgematrixlocation++]=dof;
}
nlocalNonDirEdges+=dof*dof;
}
}
}
m_comm = expList->GetComm()->GetRowComm();
m_comm->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
//Allocate arrays for block to universal map (number of expansions * p^2)
Array<OneD, long> EdgeBlockToUniversalMap(ntotaledgeentries,-1);
Array<OneD, int> localEdgeToGlobalMatrixMap(nlocalNonDirEdges,-1);
//Allocate arrays to store matrices (number of expansions * p^2)
Array<OneD, NekDouble> EdgeBlockArray(nlocalNonDirEdges,-1);
int edgematrixoffset=0;
int vGlobal;
for(n=0; n < n_exp; ++n)
{
nel = expList->GetOffset_Elmt_Id(n);
locExpansion = expList->GetExp(nel);
//loop over the edges of the expansion
for(j = 0; j < locExpansion->GetNedges(); ++j)
{
//get mesh edge id
meshEdgeId = locExpansion->as<LocalRegions::Expansion2D>()
->GetGeom2D()->GetEid(j);
nedgemodes=locExpansion->GetEdgeNcoeffs(j)-2;
if(edgeDirMap.count(meshEdgeId)==0)
{
for(k=0; k<nedgemodes*nedgemodes; ++k)
{
vGlobal=edgeglobaloffset[uniqueEdgeMap[meshEdgeId]]+k;
localEdgeToGlobalMatrixMap[edgematrixoffset+k]=vGlobal;
EdgeBlockToUniversalMap[vGlobal]
= meshEdgeId * maxEdgeDof * maxEdgeDof + k + 1;
}
edgematrixoffset+=nedgemodes*nedgemodes;
}
}
}
edgematrixoffset=0;
m_blkMat = MemoryManager<DNekBlkMat>
::AllocateSharedPtr(n_blks, n_blks, blkmatStorage);
//Here we loop over the expansion and build the block low energy
//preconditioner as well as the block versions of the transformation
//matrices.
for(cnt=n=0; n < n_exp; ++n)
{
nel = expList->GetOffset_Elmt_Id(n);
locExpansion = expList->GetExp(nel);
nVerts=locExpansion->GetGeom()->GetNumVerts();
nEdges=locExpansion->GetGeom()->GetNumEdges();
//Get statically condensed matrix
loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
//Extract boundary block (elemental S1)
bnd_mat=loc_mat->GetBlock(0,0);
//offset by number of rows
offset = bnd_mat->GetRows();
DNekScalMat &S=(*bnd_mat);
//loop over vertices of the element and return the vertex map
//for each vertex
for (v=0; v<nVerts; ++v)
{
vMap1=locExpansion->GetVertexMap(v);
//Get vertex map
globalrow = m_locToGloMap->
GetLocalToGlobalBndMap(cnt+vMap1)-nDirBnd;
if(globalrow >= 0)
{
for (m=0; m<nVerts; ++m)
{
vMap2=locExpansion->GetVertexMap(m);
//global matrix location (without offset due to
//dirichlet values)
globalcol = m_locToGloMap->
GetLocalToGlobalBndMap(cnt+vMap2)-nDirBnd;
//offset for dirichlet conditions
if (globalcol == globalrow)
{
meshVertId = locExpansion->as<LocalRegions::Expansion2D>()->GetGeom2D()->GetVid(v);
//modal connectivity between elements
sign1 = m_locToGloMap->
GetLocalToGlobalBndSign(cnt + vMap1);
sign2 = m_locToGloMap->
GetLocalToGlobalBndSign(cnt + vMap2);
vertArray[globalrow]
+= sign1*sign2*S(vMap1,vMap2);
VertBlockToUniversalMap[globalrow]
= meshVertId * maxEdgeDof * maxEdgeDof + 1;
}
}
}
}
//loop over edges of the element and return the edge map
for (eid=0; eid<nEdges; ++eid)
{
nedgemodes=locExpansion->GetEdgeNcoeffs(eid)-2;
DNekMatSharedPtr locMat =
MemoryManager<DNekMat>::AllocateSharedPtr
(nedgemodes,nedgemodes,zero,storage);
meshEdgeId = locExpansion->as<LocalRegions::Expansion2D>()->GetGeom2D()->GetEid(eid);
Array<OneD, unsigned int> edgemodearray =
locExpansion->GetEdgeInverseBoundaryMap(eid);
if(edgeDirMap.count(meshEdgeId)==0)
{
for (v=0; v<nedgemodes; ++v)
{
eMap1=edgemodearray[v];
for (m=0; m<nedgemodes; ++m)
{
eMap2=edgemodearray[m];
//modal connectivity between elements
sign1 = m_locToGloMap->
GetLocalToGlobalBndSign(cnt + eMap1);
sign2 = m_locToGloMap->
GetLocalToGlobalBndSign(cnt + eMap2);
NekDouble globalEdgeValue =
sign1*sign2*S(eMap1,eMap2);
EdgeBlockArray[edgematrixoffset+v*nedgemodes+m]=
globalEdgeValue;
}
}
edgematrixoffset+=nedgemodes*nedgemodes;
}
}
//offset for the expansion
cnt+=offset;
}
//Assemble edge matrices of each process
Array<OneD, NekDouble> GlobalEdgeBlock(ntotaledgeentries);
Vmath::Zero(ntotaledgeentries, GlobalEdgeBlock.get(), 1);
Vmath::Assmb(EdgeBlockArray.num_elements(),
EdgeBlockArray.get(),
localEdgeToGlobalMatrixMap.get(),
GlobalEdgeBlock.get());
//Exchange vertex data over different processes
if(nNonDirVerts!=0)
{
Gs::gs_data *tmp = Gs::Init(VertBlockToUniversalMap, m_comm);
Gs::Gather(vertArray, Gs::gs_add, tmp);
}
//Exchange edge data over different processes
Gs::gs_data *tmp1 = Gs::Init(EdgeBlockToUniversalMap, m_comm);
Gs::Gather(GlobalEdgeBlock, Gs::gs_add, tmp1);
// Populate vertex block
for (int i = 0; i < nNonDirVerts; ++i)
{
VertBlk->SetValue(i,i,1.0/vertArray[i]);
}
//Set the first block to be the diagonal of the vertex space
m_blkMat->SetBlock(0,0, VertBlk);
offset=0;
//Build the edge matrices from the vector
for(int loc=0; loc<nNonDirEdgeIDs; ++loc)
{
DNekMatSharedPtr gmat =
MemoryManager<DNekMat>::AllocateSharedPtr
(nedgemodes,nedgemodes,zero,storage);
for (v=0; v<nedgemodes; ++v)
{
for (m=0; m<nedgemodes; ++m)
{
NekDouble EdgeValue = GlobalEdgeBlock[offset+v*nedgemodes+m];
gmat->SetValue(v,m,EdgeValue);
}
}
m_blkMat->SetBlock(1+loc,1+loc, gmat);
offset+=edgemodeoffset[loc];
}
int totblks=m_blkMat->GetNumberOfBlockRows();
for (i=1; i< totblks; ++i)
{
unsigned int nmodes=m_blkMat->GetNumberOfRowsInBlockRow(i);
DNekMatSharedPtr tmp_mat =
MemoryManager<DNekMat>::AllocateSharedPtr
(nmodes,nmodes,zero,storage);
tmp_mat=m_blkMat->GetBlock(i,i);
tmp_mat->Invert();
m_blkMat->SetBlock(i,i,tmp_mat);
}
}
void Nektar::MultiRegions::PreconditionerBlock::BlockPreconditioner3D ( void  )
private

Definition at line 470 of file PreconditionerBlock.cpp.

References Nektar::StdRegions::StdExpansion::as(), Vmath::Assmb(), Nektar::eDIAGONAL, Nektar::SpatialDomains::eDirichlet, Nektar::eFULL, Gs::Gather(), Nektar::StdRegions::StdExpansion::GetEdgeInverseBoundaryMap(), Nektar::StdRegions::StdExpansion::GetEdgeNcoeffs(), Nektar::StdRegions::StdExpansion::GetFaceIntNcoeffs(), Nektar::StdRegions::StdExpansion::GetFaceInverseBoundaryMap(), Gs::gs_add, Gs::Init(), m_blkMat, Nektar::MultiRegions::Preconditioner::m_comm, m_linsys, m_locToGloMap, Nektar::LibUtilities::ReduceMax, and Vmath::Zero().

Referenced by v_BuildPreconditioner().

{
boost::shared_ptr<MultiRegions::ExpList>
expList=((m_linsys.lock())->GetLocMat()).lock();
LocalRegions::ExpansionSharedPtr locExpansion;
GlobalLinSysKey m_linSysKey=(m_linsys.lock())->GetKey();
StdRegions::VarCoeffMap vVarCoeffMap;
int i, j, k, nel;
int nVerts, nEdges,nFaces;
int eid, fid, n, cnt, nedgemodes, nfacemodes;
NekDouble zero = 0.0;
int vMap1, vMap2, sign1, sign2;
int m, v, eMap1, eMap2, fMap1, fMap2;
int offset, globalrow, globalcol;
//matrix storage
MatrixStorage storage = eFULL;
MatrixStorage vertstorage = eDIAGONAL;
MatrixStorage blkmatStorage = eDIAGONAL;
//local element static condensed matrices
DNekScalBlkMatSharedPtr loc_mat;
DNekScalMatSharedPtr bnd_mat;
DNekMatSharedPtr VertBlk;
int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
int nNonDirVerts = m_locToGloMap->GetNumNonDirVertexModes();
//Vertex, edge and face preconditioner matrices
VertBlk = MemoryManager<DNekMat>::
AllocateSharedPtr(nNonDirVerts,nNonDirVerts,zero,vertstorage);
Array<OneD, NekDouble> vertArray(nNonDirVerts,0.0);
Array<OneD, long> VertBlockToUniversalMap(nNonDirVerts,-1);
int n_exp = expList->GetNumElmts();
int nNonDirEdgeIDs=m_locToGloMap->GetNumNonDirEdges();
int nNonDirFaceIDs=m_locToGloMap->GetNumNonDirFaces();
//set the number of blocks in the matrix
Array<OneD,unsigned int> n_blks(1+nNonDirEdgeIDs+nNonDirFaceIDs);
n_blks[0]=nNonDirVerts;
map<int,int> edgeDirMap;
map<int,int> faceDirMap;
map<int,int> uniqueEdgeMap;
map<int,int> uniqueFaceMap;
//this should be of size total number of local edges
Array<OneD, int> edgemodeoffset(nNonDirEdgeIDs,0);
Array<OneD, int> facemodeoffset(nNonDirFaceIDs,0);
Array<OneD, int> edgeglobaloffset(nNonDirEdgeIDs,0);
Array<OneD, int> faceglobaloffset(nNonDirFaceIDs,0);
const Array<OneD, const ExpListSharedPtr>& bndCondExp = expList->GetBndCondExpansions();
StdRegions::StdExpansion2DSharedPtr bndCondFaceExp;
const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>& bndConditions = expList->GetBndConditions();
int meshVertId;
int meshEdgeId;
int meshFaceId;
const Array<OneD, const int> &extradiredges
= m_locToGloMap->GetExtraDirEdges();
for(i=0; i<extradiredges.num_elements(); ++i)
{
meshEdgeId=extradiredges[i];
edgeDirMap[meshEdgeId] = 1;
}
//Determine which boundary edges and faces have dirichlet values
for(i = 0; i < bndCondExp.num_elements(); i++)
{
cnt = 0;
for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
{
bndCondFaceExp = bndCondExp[i]->GetExp(j)->as<StdRegions::StdExpansion2D>();
if (bndConditions[i]->GetBoundaryConditionType() ==
SpatialDomains::eDirichlet)
{
for(k = 0; k < bndCondFaceExp->GetNedges(); k++)
{
meshEdgeId = bndCondFaceExp->as<LocalRegions::Expansion2D>()->GetGeom2D()->GetEid(k);
if(edgeDirMap.count(meshEdgeId) == 0)
{
edgeDirMap[meshEdgeId] = 1;
}
}
meshFaceId = bndCondFaceExp->as<LocalRegions::Expansion2D>()->GetGeom2D()->GetFid();
faceDirMap[meshFaceId] = 1;
}
}
}
int dof=0;
int maxFaceDof=0;
int maxEdgeDof=0;
int nlocalNonDirEdges=0;
int nlocalNonDirFaces=0;
int edgematrixlocation=0;
int ntotaledgeentries=0;
// Loop over all the elements in the domain and compute max edge
// DOF. Reduce across all processes to get universal maximum.
for(cnt=n=0; n < n_exp; ++n)
{
nel = expList->GetOffset_Elmt_Id(n);
locExpansion = expList->GetExp(nel);
for (j = 0; j < locExpansion->GetNedges(); ++j)
{
dof = locExpansion->GetEdgeNcoeffs(j)-2;
maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetEid(j);
if(edgeDirMap.count(meshEdgeId)==0)
{
if(uniqueEdgeMap.count(meshEdgeId)==0 && dof > 0)
{
uniqueEdgeMap[meshEdgeId]=edgematrixlocation;
edgeglobaloffset[edgematrixlocation]+=ntotaledgeentries;
edgemodeoffset[edgematrixlocation]=dof*dof;
ntotaledgeentries+=dof*dof;
n_blks[1+edgematrixlocation++]=dof;
}
nlocalNonDirEdges+=dof*dof;
}
}
}
int facematrixlocation=0;
int ntotalfaceentries=0;
// Loop over all the elements in the domain and compute max face
// DOF. Reduce across all processes to get universal maximum.
for(cnt=n=0; n < n_exp; ++n)
{
nel = expList->GetOffset_Elmt_Id(n);
locExpansion = expList->GetExp(nel);
for (j = 0; j < locExpansion->GetNfaces(); ++j)
{
dof = locExpansion->GetFaceIntNcoeffs(j);
maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetFid(j);
if(faceDirMap.count(meshFaceId)==0)
{
if(uniqueFaceMap.count(meshFaceId)==0 && dof > 0)
{
uniqueFaceMap[meshFaceId]=facematrixlocation;
facemodeoffset[facematrixlocation]=dof*dof;
faceglobaloffset[facematrixlocation]+=ntotalfaceentries;
ntotalfaceentries+=dof*dof;
n_blks[1+nNonDirEdgeIDs+facematrixlocation++]=dof;
}
nlocalNonDirFaces+=dof*dof;
}
}
}
m_comm = expList->GetComm();
m_comm->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
m_comm->AllReduce(maxFaceDof, LibUtilities::ReduceMax);
//Allocate arrays for block to universal map (number of expansions * p^2)
Array<OneD, long> EdgeBlockToUniversalMap(ntotaledgeentries,-1);
Array<OneD, long> FaceBlockToUniversalMap(ntotalfaceentries,-1);
Array<OneD, int> localEdgeToGlobalMatrixMap(nlocalNonDirEdges,-1);
Array<OneD, int> localFaceToGlobalMatrixMap(nlocalNonDirFaces,-1);
//Allocate arrays to store matrices (number of expansions * p^2)
Array<OneD, NekDouble> EdgeBlockArray(nlocalNonDirEdges,-1);
Array<OneD, NekDouble> FaceBlockArray(nlocalNonDirFaces,-1);
int edgematrixoffset=0;
int facematrixoffset=0;
int vGlobal;
for(n=0; n < n_exp; ++n)
{
nel = expList->GetOffset_Elmt_Id(n);
locExpansion = expList->GetExp(nel);
//loop over the edges of the expansion
for(j = 0; j < locExpansion->GetNedges(); ++j)
{
//get mesh edge id
meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetEid(j);
nedgemodes=locExpansion->GetEdgeNcoeffs(j)-2;
if(edgeDirMap.count(meshEdgeId)==0)
{
for(k=0; k<nedgemodes*nedgemodes; ++k)
{
vGlobal=edgeglobaloffset[uniqueEdgeMap[meshEdgeId]]+k;
localEdgeToGlobalMatrixMap[edgematrixoffset+k]=vGlobal;
EdgeBlockToUniversalMap[vGlobal]
= meshEdgeId * maxEdgeDof * maxEdgeDof + k + 1;
}
edgematrixoffset+=nedgemodes*nedgemodes;
}
}
//loop over the faces of the expansion
for(j = 0; j < locExpansion->GetNfaces(); ++j)
{
//get mesh face id
meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetFid(j);
nfacemodes = locExpansion->GetFaceIntNcoeffs(j);
//Check if face is has dirichlet values
if(faceDirMap.count(meshFaceId)==0)
{
for(k=0; k<nfacemodes*nfacemodes; ++k)
{
vGlobal=faceglobaloffset[uniqueFaceMap[meshFaceId]]+k;
localFaceToGlobalMatrixMap[facematrixoffset+k]
= vGlobal;
FaceBlockToUniversalMap[vGlobal]
= meshFaceId * maxFaceDof * maxFaceDof + k + 1;
}
facematrixoffset+=nfacemodes*nfacemodes;
}
}
}
edgematrixoffset=0;
facematrixoffset=0;
m_blkMat = MemoryManager<DNekBlkMat>
::AllocateSharedPtr(n_blks, n_blks, blkmatStorage);
//Here we loop over the expansion and build the block low energy
//preconditioner as well as the block versions of the transformation
//matrices.
for(cnt=n=0; n < n_exp; ++n)
{
nel = expList->GetOffset_Elmt_Id(n);
locExpansion = expList->GetExp(nel);
nVerts=locExpansion->GetGeom()->GetNumVerts();
nEdges=locExpansion->GetGeom()->GetNumEdges();
nFaces=locExpansion->GetGeom()->GetNumFaces();
//Get statically condensed matrix
loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
//Extract boundary block (elemental S1)
bnd_mat=loc_mat->GetBlock(0,0);
//offset by number of rows
offset = bnd_mat->GetRows();
DNekScalMat &S=(*bnd_mat);
//loop over vertices of the element and return the vertex map
//for each vertex
for (v=0; v<nVerts; ++v)
{
vMap1=locExpansion->GetVertexMap(v);
//Get vertex map
globalrow = m_locToGloMap->
GetLocalToGlobalBndMap(cnt+vMap1)-nDirBnd;
if(globalrow >= 0)
{
for (m=0; m<nVerts; ++m)
{
vMap2=locExpansion->GetVertexMap(m);
//global matrix location (without offset due to
//dirichlet values)
globalcol = m_locToGloMap->
GetLocalToGlobalBndMap(cnt+vMap2)-nDirBnd;
//offset for dirichlet conditions
if (globalcol == globalrow)
{
meshVertId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetVid(v);
//modal connectivity between elements
sign1 = m_locToGloMap->
GetLocalToGlobalBndSign(cnt + vMap1);
sign2 = m_locToGloMap->
GetLocalToGlobalBndSign(cnt + vMap2);
vertArray[globalrow]
+= sign1*sign2*S(vMap1,vMap2);
VertBlockToUniversalMap[globalrow]
= meshVertId * maxEdgeDof * maxEdgeDof + 1;
}
}
}
}
//loop over edges of the element and return the edge map
for (eid=0; eid<nEdges; ++eid)
{
nedgemodes=locExpansion->GetEdgeNcoeffs(eid)-2;
DNekMatSharedPtr locMat =
MemoryManager<DNekMat>::AllocateSharedPtr
(nedgemodes,nedgemodes,zero,storage);
meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetEid(eid);
Array<OneD, unsigned int> edgemodearray = locExpansion->GetEdgeInverseBoundaryMap(eid);
if(edgeDirMap.count(meshEdgeId)==0)
{
for (v=0; v<nedgemodes; ++v)
{
eMap1=edgemodearray[v];
for (m=0; m<nedgemodes; ++m)
{
eMap2=edgemodearray[m];
//modal connectivity between elements
sign1 = m_locToGloMap->
GetLocalToGlobalBndSign(cnt + eMap1);
sign2 = m_locToGloMap->
GetLocalToGlobalBndSign(cnt + eMap2);
NekDouble globalEdgeValue = sign1*sign2*S(eMap1,eMap2);
EdgeBlockArray[edgematrixoffset+v*nedgemodes+m]=globalEdgeValue;
}
}
edgematrixoffset+=nedgemodes*nedgemodes;
}
}
//loop over faces of the element and return the face map
for (fid=0; fid<nFaces; ++fid)
{
nfacemodes = locExpansion->GetFaceIntNcoeffs(fid);
DNekMatSharedPtr locMat =
MemoryManager<DNekMat>::AllocateSharedPtr
(nfacemodes,nfacemodes,zero,storage);
meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetFid(fid);
Array<OneD, unsigned int> facemodearray = locExpansion->GetFaceInverseBoundaryMap(fid);
if(faceDirMap.count(meshFaceId)==0)
{
for (v=0; v<nfacemodes; ++v)
{
fMap1=facemodearray[v];
for (m=0; m<nfacemodes; ++m)
{
fMap2=facemodearray[m];
//modal connectivity between elements
sign1 = m_locToGloMap->
GetLocalToGlobalBndSign(cnt + fMap1);
sign2 = m_locToGloMap->
GetLocalToGlobalBndSign(cnt + fMap2);
// Get the face-face value from the low energy matrix (S2)
NekDouble globalFaceValue = sign1*sign2*S(fMap1,fMap2);
//local face value to global face value
FaceBlockArray[facematrixoffset+v*nfacemodes+m]=globalFaceValue;
}
}
facematrixoffset+=nfacemodes*nfacemodes;
}
}
//offset for the expansion
cnt+=offset;
}
//Assemble edge matrices of each process
Array<OneD, NekDouble> GlobalEdgeBlock(ntotaledgeentries);
Vmath::Zero(ntotaledgeentries, GlobalEdgeBlock.get(), 1);
Vmath::Assmb(EdgeBlockArray.num_elements(),
EdgeBlockArray.get(),
localEdgeToGlobalMatrixMap.get(),
GlobalEdgeBlock.get());
//Assemble face matrices of each process
Array<OneD, NekDouble> GlobalFaceBlock(ntotalfaceentries);
Vmath::Zero(ntotalfaceentries, GlobalFaceBlock.get(), 1);
Vmath::Assmb(FaceBlockArray.num_elements(),
FaceBlockArray.get(),
localFaceToGlobalMatrixMap.get(),
GlobalFaceBlock.get());
//Exchange vertex data over different processes
if(nNonDirVerts!=0)
{
Gs::gs_data *tmp = Gs::Init(VertBlockToUniversalMap, m_comm);
Gs::Gather(vertArray, Gs::gs_add, tmp);
}
//Exchange edge data over different processes
Gs::gs_data *tmp1 = Gs::Init(EdgeBlockToUniversalMap, m_comm);
Gs::Gather(GlobalEdgeBlock, Gs::gs_add, tmp1);
//Exchange face data over different processes
Gs::gs_data *tmp2 = Gs::Init(FaceBlockToUniversalMap, m_comm);
Gs::Gather(GlobalFaceBlock, Gs::gs_add, tmp2);
// Populate vertex block
for (int i = 0; i < nNonDirVerts; ++i)
{
VertBlk->SetValue(i,i,1.0/vertArray[i]);
}
//Set the first block to be the diagonal of the vertex space
m_blkMat->SetBlock(0,0, VertBlk);
offset=0;
//Build the edge matrices from the vector
for(int loc=0; loc<nNonDirEdgeIDs; ++loc)
{
DNekMatSharedPtr gmat =
MemoryManager<DNekMat>::AllocateSharedPtr
(nedgemodes,nedgemodes,zero,storage);
for (v=0; v<nedgemodes; ++v)
{
for (m=0; m<nedgemodes; ++m)
{
NekDouble EdgeValue = GlobalEdgeBlock[offset+v*nedgemodes+m];
gmat->SetValue(v,m,EdgeValue);
}
}
m_blkMat->SetBlock(1+loc,1+loc, gmat);
offset+=edgemodeoffset[loc];
}
offset=0;
//Build the face matrices from the vector
for(int loc=0; loc<nNonDirFaceIDs; ++loc)
{
nfacemodes=n_blks[1+nNonDirEdgeIDs+loc];
DNekMatSharedPtr gmat =
MemoryManager<DNekMat>::AllocateSharedPtr
(nfacemodes,nfacemodes,zero,storage);
for (v=0; v<nfacemodes; ++v)
{
for (m=0; m<nfacemodes; ++m)
{
NekDouble FaceValue = GlobalFaceBlock[offset+v*nfacemodes+m];
gmat->SetValue(v,m,FaceValue);
}
}
m_blkMat->SetBlock(1+nNonDirEdgeIDs+loc,1+nNonDirEdgeIDs+loc, gmat);
offset+=facemodeoffset[loc];
}
int totblks=m_blkMat->GetNumberOfBlockRows();
for (i=1; i< totblks; ++i)
{
unsigned int nmodes=m_blkMat->GetNumberOfRowsInBlockRow(i);
DNekMatSharedPtr tmp_mat =
MemoryManager<DNekMat>::AllocateSharedPtr
(nmodes,nmodes,zero,storage);
tmp_mat=m_blkMat->GetBlock(i,i);
tmp_mat->Invert();
m_blkMat->SetBlock(i,i,tmp_mat);
}
}
void Nektar::MultiRegions::PreconditionerBlock::BlockPreconditionerHDG ( void  )
private

Definition at line 978 of file PreconditionerBlock.cpp.

References Nektar::eDIAGONAL, Gs::Gather(), Gs::gs_add, Gs::Init(), m_blkMat, m_linsys, m_locToGloMap, and Nektar::LibUtilities::ReduceMax.

Referenced by v_BuildPreconditioner().

{
boost::shared_ptr<MultiRegions::ExpList>
expList=((m_linsys.lock())->GetLocMat()).lock();
boost::shared_ptr<MultiRegions::ExpList> trace = expList->GetTrace();
LocalRegions::ExpansionSharedPtr locExpansion;
DNekScalBlkMatSharedPtr loc_mat;
DNekScalMatSharedPtr bnd_mat;
AssemblyMapDGSharedPtr asmMap = boost::dynamic_pointer_cast<
AssemblyMapDG>(m_locToGloMap);
int i, j, k, n, cnt, cnt2;
// Figure out number of Dirichlet trace elements
int nTrace = expList->GetTrace()->GetExpSize();
int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
for (cnt = n = 0; n < nTrace; ++n)
{
if (cnt >= nDir)
{
break;
}
cnt += trace->GetExp(n)->GetNcoeffs();
}
nDir = n;
// Allocate storage for block matrix. Need number of unique faces in
// trace space.
int maxTraceSize = -1;
map<int, int> arrayOffsets;
for (cnt = 0, n = nDir; n < nTrace; ++n)
{
int nCoeffs = trace->GetExp(n)->GetNcoeffs();
int nCoeffs2 = nCoeffs * nCoeffs;
arrayOffsets[n] = cnt;
cnt += nCoeffs2;
if (maxTraceSize < nCoeffs2)
{
maxTraceSize = nCoeffs2;
}
}
// Find maximum trace size.
LibUtilities::CommSharedPtr comm =
expList->GetSession()->GetComm()->GetRowComm();
comm->AllReduce(maxTraceSize, LibUtilities::ReduceMax);
// Zero matrix storage.
Array<OneD, NekDouble> tmpStore(cnt, 0.0);
// Assemble block matrices for each trace element.
for (cnt = n = 0; n < expList->GetExpSize(); ++n)
{
int elmt = expList->GetOffset_Elmt_Id(n);
locExpansion = expList->GetExp(elmt);
Array<OneD, LocalRegions::ExpansionSharedPtr> &elmtToTraceMap =
asmMap->GetElmtToTrace()[elmt];
// Block matrix (lambda)
loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
bnd_mat = loc_mat->GetBlock(0,0);
int nFacets = locExpansion->GetNumBases() == 2 ?
locExpansion->GetNedges() : locExpansion->GetNfaces();
for (cnt2 = i = 0; i < nFacets; ++i)
{
int nCoeffs = elmtToTraceMap[i]->GetNcoeffs();
int elmtId = elmtToTraceMap[i]->GetElmtId ();
// Ignore Dirichlet edges.
if (elmtId < nDir)
{
cnt += nCoeffs;
cnt2 += nCoeffs;
continue;
}
NekDouble *off = &tmpStore[arrayOffsets[elmtId]];
for (j = 0; j < nCoeffs; ++j)
{
NekDouble sign1 = asmMap->GetLocalToGlobalBndSign(
cnt + j);
for (k = 0; k < nCoeffs; ++k)
{
NekDouble sign2 = asmMap->GetLocalToGlobalBndSign(
cnt + k);
off[j*nCoeffs + k] +=
(*bnd_mat)(cnt2+j, cnt2+k) * sign1 * sign2;
}
}
cnt += nCoeffs;
cnt2 += nCoeffs;
}
}
// Set up IDs for universal numbering.
Array<OneD, long> uniIds(tmpStore.num_elements());
for (cnt = 0, n = nDir; n < nTrace; ++n)
{
LocalRegions::ExpansionSharedPtr traceExp = trace->GetExp(n);
int nCoeffs = traceExp->GetNcoeffs();
int geomId = traceExp->GetGeom()->GetGlobalID();
for (i = 0; i < nCoeffs*nCoeffs; ++i)
{
uniIds[cnt++] = geomId * maxTraceSize + i + 1;
}
}
// Assemble matrices across partitions.
Gs::gs_data *gsh = Gs::Init(uniIds, comm);
Gs::Gather(tmpStore, Gs::gs_add, gsh);
// Set up diagonal block matrix
Array<OneD, unsigned int> n_blks(nTrace - nDir);
for (n = 0; n < nTrace - nDir; ++n)
{
n_blks[n] = trace->GetExp(n + nDir)->GetNcoeffs();
}
m_blkMat = MemoryManager<DNekBlkMat>
::AllocateSharedPtr(n_blks, n_blks, eDIAGONAL);
for (n = 0; n < nTrace - nDir; ++n)
{
int nCoeffs = trace->GetExp(n + nDir)->GetNcoeffs();
DNekMatSharedPtr tmp = MemoryManager<DNekMat>
::AllocateSharedPtr(nCoeffs, nCoeffs);
NekDouble *off = &tmpStore[arrayOffsets[n+nDir]];
for (i = 0; i < nCoeffs; ++i)
{
for (j = 0; j < nCoeffs; ++j)
{
(*tmp)(i,j) = off[i*nCoeffs + j];
}
}
tmp->Invert();
m_blkMat->SetBlock(n, n, tmp);
}
}
static PreconditionerSharedPtr Nektar::MultiRegions::PreconditionerBlock::create ( const boost::shared_ptr< GlobalLinSys > &  plinsys,
const boost::shared_ptr< AssemblyMap > &  pLocToGloMap 
)
inlinestatic

Creates an instance of this class.

Definition at line 57 of file PreconditionerBlock.h.

{
PreconditionerSharedPtr p = MemoryManager<PreconditionerBlock>::AllocateSharedPtr(plinsys,pLocToGloMap);
p->InitObject();
return p;
}
void Nektar::MultiRegions::PreconditionerBlock::v_BuildPreconditioner ( )
privatevirtual

Definition at line 82 of file PreconditionerBlock.cpp.

References ASSERTL0, BlockPreconditioner2D(), BlockPreconditioner3D(), BlockPreconditionerHDG(), Nektar::StdRegions::eHybridDGHelmBndLam, Nektar::MultiRegions::GlobalMatrixKey::GetMatrixType(), and m_linsys.

{
// Different setup for HDG
GlobalLinSysKey key = m_linsys.lock()->GetKey();
if (key.GetMatrixType() == StdRegions::eHybridDGHelmBndLam)
{
BlockPreconditionerHDG();
return;
}
boost::shared_ptr<MultiRegions::ExpList>
expList=((m_linsys.lock())->GetLocMat()).lock();
StdRegions::StdExpansionSharedPtr locExpansion;
locExpansion = expList->GetExp(0);
int nDim = locExpansion->GetShapeDimension();
if (nDim == 1)
{
ASSERTL0(0,"Unknown preconditioner");
}
else if (nDim == 2)
{
BlockPreconditioner2D();
}
else if (nDim == 3)
{
BlockPreconditioner3D();
}
}
void Nektar::MultiRegions::PreconditionerBlock::v_DoPreconditioner ( const Array< OneD, NekDouble > &  pInput,
Array< OneD, NekDouble > &  pOutput 
)
privatevirtual

Definition at line 1135 of file PreconditionerBlock.cpp.

References Nektar::eWrapper, and m_locToGloMap.

{
int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
int nGlobal = m_locToGloMap->GetNumGlobalBndCoeffs();
int nNonDir = nGlobal-nDir;
DNekBlkMat &M = (*m_blkMat);
NekVector<NekDouble> r(nNonDir,pInput,eWrapper);
NekVector<NekDouble> z(nNonDir,pOutput,eWrapper);
z = M * r;
}
void Nektar::MultiRegions::PreconditionerBlock::v_InitObject ( )
privatevirtual

Definition at line 75 of file PreconditionerBlock.cpp.

References ASSERTL0, Nektar::MultiRegions::eIterativeStaticCond, and m_locToGloMap.

{
GlobalSysSolnType solvertype=m_locToGloMap->GetGlobalSysSolnType();
ASSERTL0(solvertype == MultiRegions::eIterativeStaticCond,"Solver type not valid");
}

Member Data Documentation

string Nektar::MultiRegions::PreconditionerBlock::className
static
Initial value:
GetPreconFactory().RegisterCreatorFunction(
"Block",
PreconditionerBlock::create,
"Block Diagonal Preconditioning")

Name of class.

Registers the class with the Factory.

Definition at line 68 of file PreconditionerBlock.h.

DNekBlkMatSharedPtr Nektar::MultiRegions::PreconditionerBlock::m_blkMat
protected

Definition at line 83 of file PreconditionerBlock.h.

Referenced by BlockPreconditioner2D(), BlockPreconditioner3D(), and BlockPreconditionerHDG().

const boost::weak_ptr<GlobalLinSys> Nektar::MultiRegions::PreconditionerBlock::m_linsys
protected

Definition at line 79 of file PreconditionerBlock.h.

Referenced by BlockPreconditioner2D(), BlockPreconditioner3D(), BlockPreconditionerHDG(), and v_BuildPreconditioner().

boost::shared_ptr<AssemblyMap> Nektar::MultiRegions::PreconditionerBlock::m_locToGloMap
protected

Definition at line 85 of file PreconditionerBlock.h.

Referenced by BlockPreconditioner2D(), BlockPreconditioner3D(), BlockPreconditionerHDG(), v_DoPreconditioner(), and v_InitObject().

PreconditionerType Nektar::MultiRegions::PreconditionerBlock::m_preconType
protected

Definition at line 81 of file PreconditionerBlock.h.

Generated on Wed Mar 4 2015 16:33:43 for Nektar++ by   doxygen 1.8.1.2