Commit cd6a5d31 authored by David Moxey's avatar David Moxey
Browse files

More tidying, move things up to superclass

parent 4e69d686
......@@ -74,36 +74,136 @@ namespace Nektar
const int pNumDir)
{
const int nHomDofs = pNumRows - pNumDir;
int i;
// Populate RHS vector from input
VecSetValues(m_b, nHomDofs, &m_reorderedMap[0], &pInput[pNumDir], INSERT_VALUES);
// Assemble RHS vector
VecAssemblyBegin(m_b);
VecAssemblyEnd (m_b);
// Do system solve
PetscErrorCode ierr = KSPSolve(m_ksp, m_b, m_x);
// Grab number of iterations taken
PetscInt its;
KSPGetIterationNumber(m_ksp,&its);
cout << "iteration = " << its << endl;
KSPGetIterationNumber(m_ksp, &its);
// Scatter results to local vector
VecScatterBegin(m_ctx, m_x, m_locVec, INSERT_VALUES, SCATTER_FORWARD);
VecScatterEnd (m_ctx, m_x, m_locVec, INSERT_VALUES, SCATTER_FORWARD);
// Copy results into output vector
PetscScalar *avec;
VecGetArray (m_locVec, &avec);
Vmath::Vcopy (nHomDofs, avec, 1, &pOutput[0], 1);
VecRestoreArray(m_locVec, &avec);
}
void GlobalLinSysPETSc::SetUpScatter()
{
const int nHomDofs = m_reorderedMap.size();
// Create local and global numbering systems for vector
IS isGlobal, isLocal;
ISCreateGeneral(PETSC_COMM_SELF, nHomDofs, &m_reorderedMap[0], PETSC_COPY_VALUES, &isGlobal);
ISCreateStride(PETSC_COMM_SELF, nHomDofs, 0, 1, &isLocal);
ISCreateGeneral(PETSC_COMM_SELF, nHomDofs, &m_reorderedMap[0],
PETSC_COPY_VALUES, &isGlobal);
ISCreateStride (PETSC_COMM_SELF, nHomDofs, 0, 1, &isLocal);
Vec locVec;
VecCreate (PETSC_COMM_SELF, &locVec);
VecSetSizes (locVec, nHomDofs, PETSC_DECIDE);
VecSetFromOptions(locVec);
// Create local vector for output
VecCreate (PETSC_COMM_SELF, &m_locVec);
VecSetSizes (m_locVec, m_reorderedMap.size(), PETSC_DECIDE);
VecSetFromOptions(m_locVec);
VecScatter ctx;
VecScatterCreate(m_x, isGlobal, locVec, isLocal, &ctx);
VecScatterBegin(ctx, m_x, locVec, INSERT_VALUES, SCATTER_FORWARD);
VecScatterEnd(ctx, m_x, locVec, INSERT_VALUES, SCATTER_FORWARD);
// Create scatter context
VecScatterCreate (m_x, isGlobal, m_locVec, isLocal, &m_ctx);
PetscScalar *avec;
VecGetArray(locVec, &avec);
Vmath::Vcopy(nHomDofs, avec, 1, &pOutput[0], 1);
VecRestoreArray(locVec, &avec);
// Clean up
ISDestroy(&isGlobal);
ISDestroy(&isLocal);
}
void GlobalLinSysPETSc::CalculateReordering(
const Array<OneD, const int> &glo2uniMap,
const Array<OneD, const int> &glo2unique,
const AssemblyMapSharedPtr &pLocToGloMap)
{
LibUtilities::CommSharedPtr vComm
= m_expList.lock()->GetSession()->GetComm();
const int nDirDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
const int nHomDofs = glo2uniMap.num_elements() - nDirDofs;
const int nProc = vComm->GetSize();
const int rank = vComm->GetRank();
int n, cnt;
m_nLocal = Vmath::Vsum(nHomDofs, glo2unique + nDirDofs, 1);
m_reorderedMap.resize(nHomDofs);
Array<OneD, int> localCounts(nProc, 0), localOffset(nProc, 0);
localCounts[rank] = nHomDofs;
vComm->AllReduce(localCounts, LibUtilities::ReduceSum);
for (n = 1; n < nProc; ++n)
{
localOffset[n] = localOffset[n-1] + localCounts[n-1];
}
int totHomDofs = Vmath::Vsum(nProc, localCounts, 1);
vector<unsigned int> allUniIds(totHomDofs, 0);
for (n = 0; n < nHomDofs; ++n)
{
int gid = n + nDirDofs;
allUniIds[n + localOffset[rank]] = glo2uniMap[gid];
}
vComm->AllReduce(allUniIds, LibUtilities::ReduceSum);
std::sort(allUniIds.begin(), allUniIds.end());
map<int,int> uniIdReorder;
for (cnt = n = 0; n < allUniIds.size(); ++n)
{
if (uniIdReorder.count(allUniIds[n]) > 0)
{
continue;
}
uniIdReorder[allUniIds[n]] = cnt++;
}
map<int,int>::iterator it;
for (it = uniIdReorder.begin(); it != uniIdReorder.end(); ++it)
{
cout << it->first << " -> " << it->second << endl;
}
for (n = 0; n < nHomDofs; ++n)
{
int gid = n + nDirDofs;
int uniId = glo2uniMap[gid];
ASSERTL0(uniIdReorder.count(uniId) > 0, "Error in ordering");
m_reorderedMap[n] = uniIdReorder[uniId];
}
}
void GlobalLinSysPETSc::SetUpMatVec()
{
// CREATE VECTORS
VecCreate (PETSC_COMM_WORLD, &m_x);
VecSetSizes (m_x, m_nLocal, PETSC_DECIDE);
VecSetFromOptions(m_x);
VecDuplicate (m_x, &m_b);
// CREATE MATRICES
MatCreate (PETSC_COMM_WORLD, &m_matrix);
MatSetType (m_matrix, MATMPIAIJ);
MatSetSizes (m_matrix, m_nLocal, m_nLocal,
PETSC_DETERMINE, PETSC_DETERMINE);
MatSetFromOptions(m_matrix);
MatSetUp (m_matrix);
}
}
}
......@@ -69,9 +69,18 @@ namespace Nektar
protected:
Mat m_matrix;
Vec m_x, m_b;
Vec m_x, m_b, m_locVec;
KSP m_ksp;
vector<int> m_reorderedMap;
VecScatter m_ctx;
int m_nLocal;
void SetUpScatter();
void SetUpMatVec();
void CalculateReordering(
const Array<OneD, const int> &glo2uniMap,
const Array<OneD, const int> &glo2unique,
const AssemblyMapSharedPtr &pLocToGloMap);
};
}
}
......
......@@ -68,150 +68,23 @@ namespace Nektar
"This routine should only be used when using a Full PETSc"
" matrix solve");
int i, j, n, cnt, gid1, gid2;
NekDouble sign1, sign2, value;
int nDofs = pLocToGloMap->GetNumGlobalCoeffs();
int nDirDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
int nHomDofs = nDofs - nDirDofs;
const int nDofs = pLocToGloMap->GetNumGlobalCoeffs();
const int nDirDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
const int nHomDofs = nDofs - nDirDofs;
// fill global matrix
int i, j, n, cnt, gid1, gid2, loc_lda;
NekDouble sign1, sign2, value;
DNekScalMatSharedPtr loc_mat;
LibUtilities::CommSharedPtr vComm
= pExp.lock()->GetSession()->GetComm();
const Array<OneD, const int>& glo2unique
= pLocToGloMap->GetGlobalToUniversalMapUnique();
int nLocal = Vmath::Vsum(nHomDofs, glo2unique + nDirDofs, 1);
int nProc = vComm->GetSize();
int rank = vComm->GetRank();
m_reorderedMap.resize(nHomDofs);
Array<OneD, int> localCounts(nProc, 0);
Array<OneD, int> localOffset(nProc, 0);
localCounts[rank] = nHomDofs;
vComm->AllReduce(localCounts, LibUtilities::ReduceSum);
for (i = 1; i < nProc; ++i)
{
localOffset[i] = localOffset[i-1] + localCounts[i-1];
}
int totHomDofs = Vmath::Vsum(nProc, localCounts, 1);
vector<unsigned int> allUniIds(totHomDofs, 0);
for (n = cnt = 0; n < m_expList.lock()->GetNumElmts(); ++n)
{
loc_mat = GetBlock(m_expList.lock()->GetOffset_Elmt_Id(n));
int loc_lda = loc_mat->GetRows();
// CALCULATE REORDERING MAPPING
CalculateReordering(pLocToGloMap->GetGlobalToUniversalMap(),
pLocToGloMap->GetGlobalToUniversalMapUnique(),
pLocToGloMap);
for (i = 0; i < loc_lda; ++i)
{
gid1 = pLocToGloMap->GetLocalToGlobalMap(cnt+i) - nDirDofs;
if (gid1 >= 0)
{
allUniIds[gid1+localOffset[rank]] =
pLocToGloMap->GetGlobalToUniversalMap(gid1 + nDirDofs);
}
}
cnt += loc_lda;
}
vComm->AllReduce(allUniIds, LibUtilities::ReduceSum);
std::sort(allUniIds.begin(), allUniIds.end());
map<int,int> uniIdReorder;
for (cnt = n = 0; n < allUniIds.size(); ++n)
{
if (uniIdReorder.count(allUniIds[n]) > 0)
{
continue;
}
uniIdReorder[allUniIds[n]] = cnt++;
}
for (n = cnt = 0; n < m_expList.lock()->GetNumElmts(); ++n)
{
loc_mat = GetBlock(m_expList.lock()->GetOffset_Elmt_Id(n));
int loc_lda = loc_mat->GetRows();
for (i = 0; i < loc_lda; ++i)
{
gid1 = pLocToGloMap->GetLocalToGlobalMap(cnt+i) - nDirDofs;
if (gid1 >= 0)
{
int uniid = pLocToGloMap->GetGlobalToUniversalMap(gid1 + nDirDofs);
ASSERTL0(uniIdReorder.count(uniid) > 0, "Error in ordering");
m_reorderedMap[gid1] = uniIdReorder[uniid];
}
}
cnt += loc_lda;
}
// CREATE VECTORS
VecCreate (PETSC_COMM_WORLD, &m_x);
VecSetSizes (m_x, nLocal, PETSC_DECIDE);
VecSetFromOptions(m_x);
VecDuplicate (m_x, &m_b);
// CREATE MATRICES
MatCreate (PETSC_COMM_WORLD, &m_matrix);
MatSetType (m_matrix, MATMPIAIJ);
MatSetSizes (m_matrix, nLocal, nLocal,
PETSC_DETERMINE, PETSC_DETERMINE);
MatSetFromOptions(m_matrix);
MatSetUp (m_matrix);
#if 0
// preallocate
int loc_lda;
Array<OneD, int> nnz(nDofs-NumDirBCs, 0);
vector<set<int> > hack(rows);
for(n = cnt = 0; n < m_expList.lock()->GetNumElmts(); ++n)
{
loc_mat = GetBlock(m_expList.lock()->GetOffset_Elmt_Id(n));
loc_lda = loc_mat->GetRows();
for(i = 0; i < loc_lda; ++i)
{
gid1 = pLocToGloMap->GetLocalToGlobalMap(cnt + i)-NumDirBCs;
if(gid1 < 0)
{
continue;
}
for(j = 0; j < loc_lda; ++j)
{
gid2 = pLocToGloMap->GetLocalToGlobalMap(cnt + j)
- NumDirBCs;
if(gid2 < 0)
{
continue;
}
if (hack[gid1].count(gid2) > 0)
{
continue;
}
hack[gid1].insert(gid2);
nnz[gid1]++;
}
}
cnt += loc_lda;
}
MatSeqAIJSetPreallocation(m_matrix, 0, nnz.get());
#endif
int loc_lda;
// SET UP VECTORS AND MATRIX
SetUpMatVec();
// POPULATE MATRIX
for(n = cnt = 0; n < m_expList.lock()->GetNumElmts(); ++n)
{
loc_mat = GetBlock(m_expList.lock()->GetOffset_Elmt_Id(n));
......@@ -246,6 +119,9 @@ namespace Nektar
MatAssemblyBegin(m_matrix,MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(m_matrix,MAT_FINAL_ASSEMBLY);
// SET UP SCATTER OBJECTS
SetUpScatter();
// CONSTRUCT KSP OBJECT
KSPCreate(PETSC_COMM_WORLD, &m_ksp);
KSPSetTolerances(
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment