Commit eeefef9e authored by Andrew Comerford's avatar Andrew Comerford
Browse files

Changed AssembleBnd to Scatr

parent 571819f4
......@@ -66,9 +66,6 @@ namespace Nektar
/// Global to universal unique map
Array<OneD, int> m_map;
/// Operator preconditioner matrix.
DNekMatSharedPtr m_preconditioner;
/// Tolerance of iterative solver.
NekDouble m_tolerance;
......@@ -138,8 +135,6 @@ namespace Nektar
const Array<OneD, NekDouble>& pInput,
Array<OneD, NekDouble>& pOutput) = 0;
virtual void v_ComputePreconditioner() = 0;
virtual void v_UniqueMap() = 0;
};
}
......
......@@ -176,128 +176,6 @@ namespace Nektar
}
/**
* Populates preconditioner with the identity to apply no
* preconditioning.
* @param pLocToGloMap Local to Global mapping.
*/
void GlobalLinSysIterativeFull::ComputeNullPreconditioner(
const boost::shared_ptr<AssemblyMap> &pLocToGloMap)
{
int nGlobal = pLocToGloMap->GetNumGlobalCoeffs();
int nDir = pLocToGloMap->GetNumGlobalDirBndCoeffs();
int nInt = nGlobal - nDir;
MatrixStorage storage = eDIAGONAL;
m_preconditioner = MemoryManager<DNekMat>::AllocateSharedPtr(nInt, nInt, storage);
DNekMat &M = (*m_preconditioner);
for (unsigned int i = 0; i < nInt; ++i)
{
M.SetValue(i,i,1.0);
}
}
/**
* Diagonal preconditioner computed by evaluating the local matrix
* acting on each basis vector (0,...,0,1,0,...,0). (deprecated)
* @param pLocToGloMap Local to global mapping.
*/
void GlobalLinSysIterativeFull::ComputeDiagonalPreconditioner(
const boost::shared_ptr<AssemblyMap> &pLocToGloMap)
{
boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
int nGlobal = pLocToGloMap->GetNumGlobalCoeffs();
int nLocal = pLocToGloMap->GetNumLocalCoeffs();
int nDir = pLocToGloMap->GetNumGlobalDirBndCoeffs();
int nInt = nGlobal - nDir;
MatrixStorage storage = eDIAGONAL;
m_preconditioner = MemoryManager<DNekMat>::AllocateSharedPtr(nInt, nInt, storage);
DNekMat &M = (*m_preconditioner);
Array<OneD, int> vMap = pLocToGloMap->GetLocalToGlobalMap();
for (unsigned int i = 0; i < nInt; ++i)
{
Array<OneD, NekDouble> test(nGlobal, 0.0);
Array<OneD, NekDouble> test_local(nLocal, 0.0);
test[i+nDir] = 1.0;
expList->GeneralMatrixOp(m_linSysKey, test, test, eGlobal);
M.SetValue(i,i,1.0/test[i+nDir]);
}
}
/**
* Diagonal preconditioner computed by summing the relevant elements of
* the local matrix system.
* @param pLocToGloMap Local to global mapping.
*/
void GlobalLinSysIterativeFull::ComputeDiagonalPreconditionerSum(
const boost::shared_ptr<AssemblyMap> &pLocToGloMap)
{
boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
int i,j,n,cnt,gid1,gid2;
NekDouble sign1,sign2,value;
int nGlobal = pLocToGloMap->GetNumGlobalCoeffs();
int nDir = pLocToGloMap->GetNumGlobalDirBndCoeffs();
int nInt = nGlobal - nDir;
NekDouble zero = 0.0;
// fill global matrix
DNekScalMatSharedPtr loc_mat;
Array<OneD, NekDouble> vOutput(nGlobal,0.0);
MatrixStorage storage = eDIAGONAL;
m_preconditioner = MemoryManager<DNekMat>::AllocateSharedPtr(nInt, nInt, storage);
DNekMat &M = (*m_preconditioner);
int loc_lda;
for(n = cnt = 0; n < expList->GetNumElmts(); ++n)
{
//loc_mat = vBlkMat->GetBlock(n,n);
loc_mat = GetBlock(expList->GetOffset_Elmt_Id(n));
loc_lda = loc_mat->GetRows();
for(i = 0; i < loc_lda; ++i)
{
gid1 = pLocToGloMap->GetLocalToGlobalMap(cnt + i) - nDir;
sign1 = pLocToGloMap->GetLocalToGlobalSign(cnt + i);
if(gid1 >= 0)
{
for(j = 0; j < loc_lda; ++j)
{
gid2 = pLocToGloMap->GetLocalToGlobalMap(cnt + j)
- nDir;
sign2 = pLocToGloMap->GetLocalToGlobalSign(cnt + j);
if(gid2 == gid1)
{
// When global matrix is symmetric,
// only add the value for the upper
// triangular part in order to avoid
// entries to be entered twice
value = vOutput[gid1 + nDir]
+ sign1*sign2*(*loc_mat)(i,j);
vOutput[gid1 + nDir] = value;
}
}
}
}
cnt += loc_lda;
}
// Assemble diagonal contributions across processes
pLocToGloMap->UniversalAssemble(vOutput);
// Populate preconditioner with reciprocal of diagonal elements
for (unsigned int i = 0; i < nInt; ++i)
{
M.SetValue(i,i,1.0/vOutput[i + nDir]);
}
}
/**
*
*/
......@@ -364,19 +242,6 @@ namespace Nektar
}
/**
*
*/
void GlobalLinSysIterativeFull::v_ComputePreconditioner()
{
ASSERTL1(!m_preconditioner.get(),
"Preconditioner has already been defined.");
ComputeDiagonalPreconditionerSum(m_locToGloMap);
m_map = m_locToGloMap->GetGlobalToUniversalMapUnique();
}
/**
*
*/
......
......@@ -88,26 +88,10 @@ namespace Nektar
const Array<OneD, const NekDouble> &dirForcing
= NullNekDouble1DArray);
// Populates the preconditioner with the identity matrix.
void ComputeNullPreconditioner(
const boost::shared_ptr<AssemblyMap> &pLocToGloMap);
// Populates the preconditioner with the diagonal of the operator
// matrix by applying the operator to the standard basis.
void ComputeDiagonalPreconditioner(
const boost::shared_ptr<AssemblyMap> &pLocToGloMap);
// Populates the preconditioner with the diagonal of the operator
// matrix using a more efficient summation of local matrix entries.
void ComputeDiagonalPreconditionerSum(
const boost::shared_ptr<AssemblyMap> &pLocToGloMap);
virtual void v_DoMatrixMultiply(
const Array<OneD, NekDouble>& pInput,
Array<OneD, NekDouble>& pOutput);
virtual void v_ComputePreconditioner();
virtual void v_UniqueMap();
};
......
......@@ -280,6 +280,7 @@ namespace Nektar
pLocToGloMap->GlobalToLocalBnd(F_GlobBnd,F_LocBnd);
F_LocBnd=R*F_LocBnd;
pLocToGloMap->AssembleBnd(F_LocBnd,F_HomBnd, nDirBndDofs);
//pLocToGloMap->LocalBndToGlobal(F_LocBnd,F_HomBnd, nDirBndDofs);
}
......@@ -300,24 +301,11 @@ namespace Nektar
{
DNekScalBlkMat &RT = *m_RTBlk;
NekVector<NekDouble> ml(nLocBndDofs,0.0);
NekVector<NekDouble> MultVector(nGlobHomBndDofs,1.0);
m_locToGloMap->GlobalToLocalBnd(MultVector,ml,nDirBndDofs);
m_locToGloMap->AssembleBnd(ml,MultVector,nDirBndDofs);
for(int i=0; i<nGlobHomBndDofs; ++i)
{
MultVector[i]=1/MultVector[i];
}
pLocToGloMap->GlobalToLocalBnd(V_GlobHomBnd,V_LocBnd, nDirBndDofs);
V_LocBnd=RT*V_LocBnd;
pLocToGloMap->AssembleBnd(V_LocBnd,V_GlobHomBnd, nDirBndDofs);
V_GlobHomBnd=V_GlobHomBnd*MultVector;
pLocToGloMap->LocalBndToGlobal(V_LocBnd,V_GlobHomBnd, nDirBndDofs);
}
}
else
......@@ -852,11 +840,6 @@ namespace Nektar
}
}
void GlobalLinSysIterativeStaticCond::v_ComputePreconditioner()
{
//
}
void GlobalLinSysIterativeStaticCond::v_UniqueMap()
{
m_map = m_locToGloMap->GetGlobalToUniversalBndMapUnique();
......
......@@ -149,18 +149,11 @@ namespace Nektar
void ConstructNextLevelCondensedSystem(
const boost::shared_ptr<AssemblyMap>& locToGloMap);
/// Compute a diagonal preconditioner of the Shur-complement matrix.
void ComputeDiagonalPreconditioner(
const boost::shared_ptr<AssemblyMap> &pLocToGloMap);
/// Perform a Shur-complement matrix multiply operation.
virtual void v_DoMatrixMultiply(
const Array<OneD, NekDouble>& pInput,
Array<OneD, NekDouble>& pOutput);
/// Compute the preconditioner for the Shur-complement matrix.
virtual void v_ComputePreconditioner();
virtual void v_UniqueMap();
};
......
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