#include <PreconditionerBlock.h>

Inheritance diagram for Nektar::MultiRegions::PreconditionerBlock:

Collaboration diagram for Nektar::MultiRegions::PreconditionerBlock:

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. More...

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. More...

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

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	BlockPreconditionerCG (void)
	Construct a block preconditioner from $\mathbf{S}_{1}$ for the continuous Galerkin system. More...

void	BlockPreconditionerHDG (void)
	Construct a block preconditioner for the hybridized discontinuous Galerkin system. More...

virtual void	v_InitObject ()

virtual void	v_DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
	Apply preconditioner to `pInput` and store the result in `pOutput`. More...

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}$ . More...

Detailed Description

Definition at line 52 of file PreconditionerBlock.h.

Constructor & Destructor Documentation

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

Definition at line 67 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.

75 {}

Member Function Documentation

void Nektar::MultiRegions::PreconditionerBlock::BlockPreconditionerCG ( void )

private

Construct a block preconditioner from $\mathbf{S}_{1}$ for the continuous Galerkin system.

The preconditioner is defined as:

$\mathbf{M}^{-1}=\left[\begin{array}{ccc} \mathrm{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 and the subscript denotes the portion of the Schur complement associated with the vertex, edge and face blocks respectively.

Definition at line 116 of file PreconditionerBlock.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL1, Nektar::MultiRegions::DeterminePeriodicEdgeOrientId(), Nektar::MultiRegions::DeterminePeriodicFaceOrient(), Nektar::eDIAGONAL, Gs::Gather(), Gs::gs_add, Gs::Init(), Nektar::iterator, m_blkMat, Nektar::MultiRegions::Preconditioner::m_comm, m_linsys, m_locToGloMap, Nektar::LibUtilities::ReduceMax, sign, and Vmath::Vadd().

Referenced by v_BuildPreconditioner().

         {
             ExpListSharedPtr expList = m_linsys.lock()->GetLocMat().lock();
             LocalRegions::ExpansionSharedPtr exp;
             DNekScalBlkMatSharedPtr loc_mat;
             DNekScalMatSharedPtr    bnd_mat;
 
             int nel, i, j, k, n, cnt, gId;
             int meshVertId, meshEdgeId, meshFaceId;
 
             const int nExp = expList->GetExpSize();
             const int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
 
             // Grab periodic geometry information.
             PeriodicMap::iterator pIt;
             PeriodicMap periodicVerts, periodicEdges, periodicFaces;
             expList->GetPeriodicEntities(
                 periodicVerts, periodicEdges, periodicFaces);
 
             // The vectors below are of size 3 to have separate storage for
             // vertices, edges and faces.
 
             // Maps from geometry ID to the matrix representing the extracted
             // portion of S_1. For example idMats[2] folds the S_1 face blocks.
             vector<map<int, vector<NekDouble> > > idMats(3);
 
             // Maps from the global ID, as obtained from AssemblyMapCG's
             // localToGlobalMap, to the geometry ID.
             vector<map<int, int> > gidMeshIds(3);
 
             // Maps from the global ID to the number of degrees of freedom for
             // this geometry object.
             vector<map<int, int> > gidDofs(3);
 
             // Array containing maximum information needed for the universal
             // numbering later. For i = 0,1,2 for each geometry dimension:
             //
             // maxVertIds[2*i]   = maximum geometry ID at dimension i
             // maxVertIds[2*i+1] = maximum number of degrees of freedom for all
             //                     elements of dimension i.
             Array<OneD, int> maxVertIds(6, -1);
 
             // Iterator for idMats members.
             map<int, vector<NekDouble> >::iterator gIt;
 
             // Figure out mapping from each elemental contribution to offset in
             // (vert,edge,face) triples.
             for (cnt = n = 0; n < nExp; ++n)
             {
                 nel = expList->GetOffset_Elmt_Id(n);
                 exp = expList->GetExp(nel);
 
                 // Grab reference to local Schur complement matrix.
                 DNekScalMatSharedPtr schurMat =
                     m_linsys.lock()->GetStaticCondBlock(n)->GetBlock(0,0);
 
                 // Process vertices to extract relevant portion of the Schur
                 // complement matrix.
                 for (i = 0; i < exp->GetNverts(); ++i)
                 {
                     meshVertId = exp->GetGeom()->GetVid(i);
                     int locId = exp->GetVertexMap(i);
 
                     // Get the global ID of this vertex.
                     gId = m_locToGloMap->GetLocalToGlobalMap(
                         cnt + locId) - nDirBnd;
 
                     // Ignore all Dirichlet vertices.
                     if (gId < 0)
                     {
                         continue;
                     }
 
                     gidDofs[0][gId] = 1;
 
                     // Extract vertex value from Schur complement matrix.
                     NekDouble vertVal = (*schurMat)(locId,locId);
 
                     // See if we have processed this vertex from another
                     // element.
                     gIt = idMats[0].find(gId);
 
                     if (gIt == idMats[0].end())
                     {
                         // If not then put our 'matrix' inside idMats.
                         idMats[0][gId] = vector<NekDouble>(1, vertVal);
                     }
                     else
                     {
                         // Otherwise combine with the value that is already
                         // there (i.e. do assembly on this degree of freedom).
                         gIt->second[0] += vertVal;
                     }
 
                     // Now check to see if the vertex is periodic. If it is,
                     // then we change meshVertId to be the minimum of all the
                     // other periodic vertices, so that we don't end up
                     // duplicating the matrix in our final block matrix.
                     pIt = periodicVerts.find(meshVertId);
                     if (pIt != periodicVerts.end())
                     {
                         for (j = 0; j < pIt->second.size(); ++j)
                         {
                             meshVertId = min(meshVertId, pIt->second[j].id);
                         }
                     }
 
                     // Finally record the other information we need into the
                     // other matrices.
                     gidMeshIds[0][gId] = meshVertId;
                     maxVertIds[0] = max(maxVertIds[0], meshVertId);
                     maxVertIds[1] = 1;
                 }
 
                 // Process edges. This logic is mostly the same as the previous
                 // block.
                 for (i = 0; i < exp->GetNedges(); ++i)
                 {
                     meshEdgeId = exp->GetGeom()->GetEid(i);
 
                     Array<OneD, unsigned int> bmap, bmap2;
                     Array<OneD, int> sign;
                     StdRegions::Orientation edgeOrient =
                         exp->GetGeom()->GetEorient(i);
 
                     // Check if this edge is periodic. We may need to flip
                     // orientation if it is.
                     pIt = periodicEdges.find(meshEdgeId);
                     if (pIt != periodicEdges.end())
                     {
                         pair<int, StdRegions::Orientation> idOrient =
                             DeterminePeriodicEdgeOrientId(
                                 meshEdgeId, edgeOrient, pIt->second);
                         meshEdgeId = idOrient.first;
                         edgeOrient = idOrient.second;
                     }
 
                     // Grab edge interior map, and the edge inverse boundary
                     // map, so that we can extract this edge from the Schur
                     // complement matrix.
                     exp->GetEdgeInteriorMap(i, edgeOrient, bmap, sign);
                     bmap2 = exp->GetEdgeInverseBoundaryMap(i);
 
                     // Allocate temporary storage for the extracted edge matrix.
                     const int nEdgeCoeffs = bmap.num_elements();
                     vector<NekDouble> tmpStore(nEdgeCoeffs*nEdgeCoeffs);
 
                     gId = m_locToGloMap->GetLocalToGlobalMap(cnt + bmap[0]);
 
                     for (j = 0; j < nEdgeCoeffs; ++j)
                     {
                         // We record the minimum ID from the edge for our
                         // maps. This follows the logic that the assembly map
                         // ordering will always give us a contiguous ordering of
                         // global degrees of freedom for edge interior
                         // coefficients.
                         gId = min(gId,
                                   m_locToGloMap->GetLocalToGlobalMap(
                                       cnt + bmap[j])
                                   - nDirBnd);
 
                         // Ignore Dirichlet edges.
                         if (gId < 0)
                         {
                             continue;
                         }
 
                         const NekDouble sign1 = sign[j];
 
                         // Extract this edge, along with sign array for assembly
                         // later.
                         for (k = 0; k < nEdgeCoeffs; ++k)
                         {
                             tmpStore[k+j*nEdgeCoeffs] =
                                 sign1*sign[k]*(*schurMat)(bmap2[j], bmap2[k]);
                         }
                     }
 
                     if (gId < 0)
                     {
                         continue;
                     }
 
                     gidDofs[1][gId] = nEdgeCoeffs;
 
                     // Assemble this edge matrix with another one, if it exists.
                     gIt = idMats[1].find(gId);
                     if (gIt == idMats[1].end())
                     {
                         idMats[1][gId] = tmpStore;
                     }
                     else
                     {
                         ASSERTL1(tmpStore.size() == gIt->second.size(),
                                  "Number of modes mismatch");
                         Vmath::Vadd(nEdgeCoeffs*nEdgeCoeffs, &gIt->second[0], 1,
                                     &tmpStore[0], 1, &gIt->second[0], 1);
                     }
 
                     gidMeshIds[1][gId] = meshEdgeId;
                     maxVertIds[2] = max(maxVertIds[2], meshEdgeId);
                     maxVertIds[3] = max(maxVertIds[3], nEdgeCoeffs);
                 }
 
                 // Process faces. This logic is mostly the same as the previous
                 // block.
                 for (i = 0; i < exp->GetNfaces(); ++i)
                 {
                     meshFaceId = exp->GetGeom()->GetFid(i);
 
                     Array<OneD, unsigned int> bmap, bmap2;
                     Array<OneD, int> sign;
                     StdRegions::Orientation faceOrient =
                         exp->GetGeom()->GetForient(i);
 
                     // Check if this face is periodic. We may need to flip
                     // orientation if it is.
                     pIt = periodicFaces.find(meshFaceId);
                     if (pIt != periodicFaces.end())
                     {
                         meshFaceId = min(meshFaceId, pIt->second[0].id);
                         faceOrient = DeterminePeriodicFaceOrient(
                             faceOrient, pIt->second[0].orient);
                     }
 
                     exp->GetFaceInteriorMap(i, faceOrient, bmap, sign);
                     bmap2 = exp->GetFaceInverseBoundaryMap(i);
 
                     // Allocate temporary storage for the extracted face matrix.
                     const int nFaceCoeffs = bmap.num_elements();
                     vector<NekDouble> tmpStore(nFaceCoeffs*nFaceCoeffs);
 
                     gId = m_locToGloMap->GetLocalToGlobalMap(cnt + bmap[0]);
 
                     for (j = 0; j < nFaceCoeffs; ++j)
                     {
                         gId = min(gId,
                                   m_locToGloMap->GetLocalToGlobalMap(
                                       cnt + bmap[j])
                                   - nDirBnd);
 
                         // Ignore Dirichlet faces.
                         if (gId < 0)
                         {
                             continue;
                         }
 
                         const NekDouble sign1 = sign[j];
 
                         // Extract this face, along with sign array for assembly
                         // later.
                         for (k = 0; k < nFaceCoeffs; ++k)
                         {
                             tmpStore[k+j*nFaceCoeffs] =
                                 sign1*sign[k]*(*schurMat)(bmap2[j], bmap2[k]);
                         }
                     }
 
                     if (gId < 0)
                     {
                         continue;
                     }
 
                     gidDofs[2][gId] = nFaceCoeffs;
 
                     // Assemble this face matrix with another one, if it exists.
                     gIt = idMats[2].find(gId);
                     if (gIt == idMats[2].end())
                     {
                         idMats[2][gId] = tmpStore;
                     }
                     else
                     {
                         ASSERTL1(tmpStore.size() == gIt->second.size(),
                                  "Number of modes mismatch");
                         Vmath::Vadd(nFaceCoeffs*nFaceCoeffs, &gIt->second[0], 1,
                                     &tmpStore[0], 1, &gIt->second[0], 1);
                     }
 
                     gidMeshIds[2][gId] = meshFaceId;
                     maxVertIds[4] = max(maxVertIds[4], meshFaceId);
                     maxVertIds[5] = max(maxVertIds[5], nFaceCoeffs);
                 }
 
                 cnt += exp->GetNcoeffs();
             }
 
             // Perform a reduction to find maximum vertex, edge and face
             // geometry IDs.
             m_comm = expList->GetSession()->GetComm()->GetRowComm();
             m_comm->AllReduce(maxVertIds, LibUtilities::ReduceMax);
 
             // Concatenate all matrices into contiguous storage and figure out
             // universal ID numbering.
             vector<NekDouble> storageBuf;
             vector<long> globalToUniversal;
 
             for (i = 0, cnt = 1; i < 3; ++i)
             {
                 const int maxDofs = maxVertIds[2*i+1];
 
                 // Note that iterating over the map uses the property that keys
                 // are accessed in order of ascending order, putting everything
                 // in the right order for the global system.
                 for (gIt = idMats[i].begin(); gIt != idMats[i].end(); ++gIt)
                 {
                     // Copy matrix into storage.
                     storageBuf.insert(storageBuf.end(),
                                       gIt->second.begin(), gIt->second.end());
 
                     // Get mesh ID from global ID number.
                     ASSERTL1(gidMeshIds[i].count(gIt->first) > 0,
                              "Unable to find global ID " +
                              boost::lexical_cast<string>(gIt->first) +
                              " inside map");
                     meshVertId = gidMeshIds[i][gIt->first];
 
                     for (j = 0; j < gIt->second.size(); ++j)
                     {
                         globalToUniversal.push_back(
                             cnt + meshVertId*maxDofs*maxDofs + j);
                     }
 
                     // Free up the temporary storage.
                     gIt->second.clear();
                 }
 
                 cnt += (maxVertIds[2*i]+1)*maxDofs*maxDofs;
             }
 
             ASSERTL1(storageBuf.size() == globalToUniversal.size(),
                      "Storage buffer and global to universal map size does "
                      "not match");
 
             Array<OneD, NekDouble> storageData(
                 storageBuf.size(), &storageBuf[0]);
             Array<OneD, long> globalToUniversalMap(
                 globalToUniversal.size(), &globalToUniversal[0]);
 
             // Use GS to assemble data between processors.
             Gs::gs_data *tmpGs = Gs::Init(
                 globalToUniversalMap, m_comm,
                 expList->GetSession()->DefinesCmdLineArgument("verbose"));
             Gs::Gather(storageData, Gs::gs_add, tmpGs);
 
             // Figure out what storage we need in the block matrix.
             Array<OneD, unsigned int> n_blks(
                 1 + idMats[1].size() + idMats[2].size());
 
             // Vertex block is a diagonal matrix.
             n_blks[0] = idMats[0].size();
 
             // Now extract number of rows in each edge and face block from the
             // gidDofs map.
             cnt = 1;
             for (i = 1; i < 3; ++i)
             {
                 for (gIt = idMats[i].begin(); gIt != idMats[i].end(); ++gIt)
                 {
                     ASSERTL1(gidDofs[i].count(gIt->first) > 0,
                              "Unable to find number of degrees of freedom for "
                              "global ID " + boost::lexical_cast<string>(
                                  gIt->first));
 
                     n_blks[cnt++] = gidDofs[i][gIt->first];
                 }
             }
 
             // Allocate storage for the block matrix.
             m_blkMat = MemoryManager<DNekBlkMat>
                 ::AllocateSharedPtr(n_blks, n_blks, eDIAGONAL);
 
             // We deal with the vertex matrix separately since all vertices can
             // be condensed into a single, block-diagonal matrix.
             DNekMatSharedPtr vertMat = MemoryManager<DNekMat>
                 ::AllocateSharedPtr(n_blks[0], n_blks[0], 0.0, eDIAGONAL);
 
             // Fill the vertex matrix with the inverse of each vertex value.
             cnt = 0;
             for (gIt = idMats[0].begin(); gIt != idMats[0].end(); ++gIt, ++cnt)
             {
                 (*vertMat)(cnt, cnt) = 1.0/storageData[cnt];
             }
 
             // Put the vertex matrix in the block matrix.
             m_blkMat->SetBlock(0,0,vertMat);
 
             // Finally, grab the matrices from the block storage, invert them
             // and place them in the correct position inside the block matrix.
             int cnt2 = 1;
             for (i = 1; i < 3; ++i)
             {
                 for (gIt = idMats[i].begin(); gIt != idMats[i].end();
                      ++gIt, ++cnt2)
                 {
                     int nDofs = gidDofs[i][gIt->first];
 
                     DNekMatSharedPtr tmp = MemoryManager<DNekMat>
                         ::AllocateSharedPtr(nDofs, nDofs);
 
                     for (j = 0; j < nDofs; ++j)
                     {
                         for (k = 0; k < nDofs; ++k)
                         {
                             (*tmp)(j,k) = storageData[k+j*nDofs + cnt];
                         }
                     }
 
                     cnt += nDofs*nDofs;
 
                     tmp->Invert();
                     m_blkMat->SetBlock(cnt2, cnt2, tmp);
                 }
             }
         }

void Nektar::MultiRegions::PreconditionerBlock::BlockPreconditionerHDG ( void )

private

Construct a block preconditioner for the hybridized discontinuous Galerkin system.

This system is built in a similar fashion to its continuous variant found in PreconditionerBlock::BlockPreconditionerCG. In this setting however, the matrix is constructed as:

$M^{-1} = \mathrm{Diag}[ (\mathbf{S_{1}})_{f}^{-1} ]$

where each matrix is the Schur complement system restricted to a single face of the trace system.

Definition at line 545 of file PreconditionerBlock.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eDIAGONAL, Gs::Gather(), Gs::gs_add, Gs::Init(), m_blkMat, Nektar::MultiRegions::Preconditioner::m_comm, 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.
             m_comm = expList->GetSession()->GetComm()->GetRowComm();
             m_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, m_comm,
                 expList->GetSession()->DefinesCmdLineArgument("verbose"));
             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 56 of file PreconditionerBlock.h.

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

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

void Nektar::MultiRegions::PreconditionerBlock::v_BuildPreconditioner ( )

privatevirtual

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 85 of file PreconditionerBlock.cpp.

References BlockPreconditionerCG(), BlockPreconditionerHDG(), Nektar::StdRegions::eHybridDGHelmBndLam, Nektar::MultiRegions::GlobalMatrixKey::GetMatrixType(), and m_linsys.

         {
             GlobalLinSysKey key = m_linsys.lock()->GetKey();
 
             // Different setup for HDG and CG.
             if (key.GetMatrixType() == StdRegions::eHybridDGHelmBndLam)
             {
                 BlockPreconditionerHDG();
             }
             else
             {
                 BlockPreconditionerCG();
             }
         }

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

privatevirtual

Apply preconditioner to pInput and store the result in pOutput.

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 704 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

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 77 of file PreconditionerBlock.cpp.

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

         {
             GlobalSysSolnType solvertype=m_locToGloMap->GetGlobalSysSolnType();
             ASSERTL0(solvertype == MultiRegions::eIterativeStaticCond ||
                      solvertype == MultiRegions::ePETScStaticCond,
                      "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 80 of file PreconditionerBlock.h.

Referenced by BlockPreconditionerCG(), and BlockPreconditionerHDG().

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

protected

Definition at line 78 of file PreconditionerBlock.h.

Referenced by BlockPreconditionerCG(), BlockPreconditionerHDG(), and v_BuildPreconditioner().

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

protected

Definition at line 81 of file PreconditionerBlock.h.

Referenced by BlockPreconditionerCG(), BlockPreconditionerHDG(), v_DoPreconditioner(), and v_InitObject().

PreconditionerType Nektar::MultiRegions::PreconditionerBlock::m_preconType

protected

Definition at line 79 of file PreconditionerBlock.h.

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Attributes

Private Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation