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

Working version with low energy precon

parent e8ac039d
......@@ -161,8 +161,39 @@ namespace Nektar
// Allocate memory for top-level structure
SetupTopLevel(m_locToGloMap);
// Setup Block Matrix systems
int n, n_exp = m_expList.lock()->GetNumElmts();
MatrixStorage blkmatStorage = eDIAGONAL;
const Array<OneD,const unsigned int>& nbdry_size
= m_locToGloMap->GetNumLocalBndCoeffsPerPatch();
m_S1Blk = MemoryManager<DNekScalBlkMat>
::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
// Preserve original matrix in m_S1Blk
for (n = 0; n < n_exp; ++n)
{
DNekScalMatSharedPtr mat = m_schurCompl->GetBlock(n, n);
m_S1Blk->SetBlock(n, n, mat);
}
// Build preconditioner
m_precon->BuildPreconditioner();
// Do transform of Schur complement matrix
for (n = 0; n < n_exp; ++n)
{
if (m_linSysKey.GetMatrixType() !=
StdRegions::eHybridDGHelmBndLam)
{
DNekScalMatSharedPtr mat = m_S1Blk->GetBlock(n, n);
DNekScalMatSharedPtr t = m_precon->TransformedSchurCompl(
m_expList.lock()->GetOffset_Elmt_Id(n), mat);
m_schurCompl->SetBlock(n, n, t);
}
}
// Construct this level
Initialise(m_locToGloMap);
}
......@@ -175,7 +206,6 @@ namespace Nektar
}
#if 0
DNekScalBlkMatSharedPtr GlobalLinSysIterativeStaticCond::
v_GetStaticCondBlock(unsigned int n)
{
......@@ -193,7 +223,7 @@ namespace Nektar
return schurComplBlock;
}
#endif
/**
* Assemble the schur complement matrix from the block matrices stored
* in #m_blkMatrices and the given local to global mapping information.
......@@ -202,6 +232,9 @@ namespace Nektar
void GlobalLinSysIterativeStaticCond::v_AssembleSchurComplement(
const AssemblyMapSharedPtr pLocToGloMap)
{
int i,j,n,cnt,gid1,gid2;
NekDouble sign1,sign2;
bool doGlobalOp = m_expList.lock()->GetGlobalOptParam()->
DoGlobalMatOp(m_linSysKey.GetMatrixType());
......@@ -213,9 +246,6 @@ namespace Nektar
return;
}
int i,j,n,cnt,gid1,gid2;
NekDouble sign1,sign2;
int nBndDofs = pLocToGloMap->GetNumGlobalBndCoeffs();
int NumDirBCs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
unsigned int rows = nBndDofs - NumDirBCs;
......@@ -433,16 +463,14 @@ namespace Nektar
{
int nLocal = m_locToGloMap->GetNumLocalBndCoeffs();
int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
bool doGlobalOp = m_expList.lock()->GetGlobalOptParam()->
DoGlobalMatOp(m_linSysKey.GetMatrixType());
if(doGlobalOp)
{
// Do matrix multiply globally
Array<OneD, NekDouble> in = pInput + nDir;
Array<OneD, NekDouble> out = pOutput+ nDir;
Array<OneD, NekDouble> in = pInput + nDir;
Array<OneD, NekDouble> out = pOutput + nDir;
m_sparseSchurCompl->Multiply(in,out);
m_locToGloMap->UniversalAssembleBnd(pOutput, nDir);
......@@ -479,6 +507,34 @@ namespace Nektar
m_map = m_locToGloMap->GetGlobalToUniversalBndMapUnique();
}
DNekScalBlkMatSharedPtr GlobalLinSysIterativeStaticCond::v_PreSolve(
int scLevel,
NekVector<NekDouble> &F_GlobBnd)
{
if (scLevel == 0)
{
Set_Rhs_Magnitude(F_GlobBnd);
return m_S1Blk;
}
else
{
return m_schurCompl;
}
}
void GlobalLinSysIterativeStaticCond::v_BasisTransform(
Array<OneD, NekDouble>& pInOut,
int offset)
{
m_precon->DoTransformToLowEnergy(pInOut, offset);
}
void GlobalLinSysIterativeStaticCond::v_BasisInvTransform(
Array<OneD, NekDouble>& pInOut)
{
m_precon->DoTransformFromLowEnergy(pInOut);
}
GlobalLinSysStaticCondSharedPtr GlobalLinSysIterativeStaticCond::v_Recurse(
const GlobalLinSysKey &mkey,
const boost::weak_ptr<ExpList> &pExpList,
......
......@@ -115,7 +115,7 @@ namespace Nektar
virtual ~GlobalLinSysIterativeStaticCond();
protected:
//virtual DNekScalBlkMatSharedPtr v_GetStaticCondBlock(unsigned int n);
virtual DNekScalBlkMatSharedPtr v_GetStaticCondBlock(unsigned int n);
virtual GlobalLinSysStaticCondSharedPtr v_Recurse(
const GlobalLinSysKey &mkey,
const boost::weak_ptr<ExpList> &pExpList,
......@@ -125,6 +125,15 @@ namespace Nektar
const DNekScalBlkMatSharedPtr pInvD,
const boost::shared_ptr<AssemblyMap> &locToGloMap);
virtual DNekScalBlkMatSharedPtr v_PreSolve(
int scLevel,
NekVector<NekDouble> &F_GlobBnd);
virtual void v_BasisTransform(
Array<OneD, NekDouble>& pInOut,
int offset);
virtual void v_BasisInvTransform(
Array<OneD, NekDouble>& pInOut);
private:
DNekScalBlkMatSharedPtr m_S1Blk;
/// Dense storage for block Schur complement matrix
......
......@@ -229,7 +229,7 @@ namespace Nektar
nGlobBndDofs, F, pert, pLocToGloMap, nDirBndDofs);
// Transform back to original basis
v_BasisInvTransform(F, nDirBndDofs);
v_BasisInvTransform(pert);
// Add back initial conditions onto difference
Vmath::Vadd(nGlobHomBndDofs,&out[nDirBndDofs],1,
......@@ -343,7 +343,6 @@ namespace Nektar
DNekScalBlkMatSharedPtr loc_schur
= GlobalLinSys::v_GetStaticCondBlock(
m_expList.lock()->GetOffset_Elmt_Id(n));
DNekScalMatSharedPtr t;
m_schurCompl->SetBlock(n, n, t = loc_schur->GetBlock(0,0));
m_BinvD ->SetBlock(n, n, t = loc_schur->GetBlock(0,1));
......
......@@ -79,8 +79,7 @@ namespace Nektar
}
virtual void v_BasisInvTransform(
Array<OneD, NekDouble>& pInOut,
int offset)
Array<OneD, NekDouble>& pInOut)
{
}
......
......@@ -178,8 +178,8 @@ namespace Nektar
* \brief Get block elemental transposed transformation matrix
* \f$\mathbf{R}^{T}\f$
*/
DNekScalBlkMatSharedPtr Preconditioner::v_TransformedSchurCompl(
int offset, const boost::shared_ptr<DNekScalBlkMat> &loc_mat)
DNekScalMatSharedPtr Preconditioner::v_TransformedSchurCompl(
int offset, const boost::shared_ptr<DNekScalMat> &loc_mat)
{
return loc_mat;
}
......
......@@ -127,8 +127,8 @@ namespace Nektar
inline const DNekScalBlkMatSharedPtr&
GetBlockTransposedTransformationMatrix() const;
inline DNekScalBlkMatSharedPtr TransformedSchurCompl(
int offset, const boost::shared_ptr<DNekScalBlkMat > &loc_mat);
inline DNekScalMatSharedPtr TransformedSchurCompl(
int offset, const boost::shared_ptr<DNekScalMat > &loc_mat);
protected:
const boost::weak_ptr<GlobalLinSys> m_linsys;
......@@ -137,8 +137,8 @@ namespace Nektar
boost::shared_ptr<AssemblyMap> m_locToGloMap;
LibUtilities::CommSharedPtr m_comm;
virtual DNekScalBlkMatSharedPtr v_TransformedSchurCompl(
int offset, const boost::shared_ptr<DNekScalBlkMat > &loc_mat);
virtual DNekScalMatSharedPtr v_TransformedSchurCompl(
int offset, const boost::shared_ptr<DNekScalMat > &loc_mat);
private:
......@@ -195,8 +195,8 @@ namespace Nektar
/**
*
*/
inline DNekScalBlkMatSharedPtr Preconditioner::TransformedSchurCompl(
int offset, const boost::shared_ptr<DNekScalBlkMat > &loc_mat)
inline DNekScalMatSharedPtr Preconditioner::TransformedSchurCompl(
int offset, const boost::shared_ptr<DNekScalMat > &loc_mat)
{
return v_TransformedSchurCompl(offset,loc_mat);
}
......
......@@ -96,10 +96,10 @@ namespace Nektar
}
DNekScalBlkMatSharedPtr PreconditionerLinearWithLowEnergy::
v_TransformedSchurCompl(int offset, const boost::shared_ptr<DNekScalBlkMat > &loc_mat)
DNekScalMatSharedPtr PreconditionerLinearWithLowEnergy::
v_TransformedSchurCompl(int offset, const boost::shared_ptr<DNekScalMat > &loc_mat)
{
DNekScalBlkMatSharedPtr returnval;
DNekScalMatSharedPtr returnval;
returnval=m_lowEnergyPrecon->TransformedSchurCompl(offset,loc_mat);
return returnval;
}
......
......@@ -87,8 +87,8 @@ namespace Nektar
virtual void v_DoTransformFromLowEnergy(
Array<OneD, NekDouble>& pInput);
virtual DNekScalBlkMatSharedPtr
v_TransformedSchurCompl(int offset, const boost::shared_ptr<DNekScalBlkMat > &loc_mat);
virtual DNekScalMatSharedPtr
v_TransformedSchurCompl(int offset, const boost::shared_ptr<DNekScalMat > &loc_mat);
virtual void v_DoPreconditioner(
const Array<OneD, NekDouble>& pInput,
......
......@@ -1212,10 +1212,10 @@ namespace Nektar
* the low energy equivalent
* i.e. \f$\mathbf{S}_{2}=\mathbf{R}\mathbf{S}_{1}\mathbf{R}^{T}\f$
*/
DNekScalBlkMatSharedPtr PreconditionerLowEnergy::
DNekScalMatSharedPtr PreconditionerLowEnergy::
v_TransformedSchurCompl(
int offset,
const boost::shared_ptr<DNekScalBlkMat > &loc_mat)
const boost::shared_ptr<DNekScalMat > &loc_mat)
{
boost::shared_ptr<MultiRegions::ExpList>
expList=((m_linsys.lock())->GetLocMat()).lock();
......@@ -1227,7 +1227,7 @@ namespace Nektar
unsigned int nint=ncoeffs-nbnd;
//This is the SC elemental matrix in the orginal basis (S1)
DNekScalMatSharedPtr pS1=loc_mat->GetBlock(0,0);
DNekScalMatSharedPtr pS1=loc_mat;
//Transformation matrices
map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transmatrixmap;
......@@ -1261,16 +1261,10 @@ namespace Nektar
RS1=R*S1;
S2=RS1*RT;
DNekScalBlkMatSharedPtr returnval;
DNekScalMatSharedPtr tmp_mat;
unsigned int exp_size[] = {nbnd, nint};
unsigned int nblks = 1;
returnval = MemoryManager<DNekScalBlkMat>::
AllocateSharedPtr(nblks, nblks, exp_size, exp_size);
tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, pS2);
returnval->SetBlock(0,0,tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0,pS2));
return returnval;
return tmp_mat;
}
/**
......
......@@ -151,8 +151,8 @@ namespace Nektar
virtual void v_BuildPreconditioner();
virtual DNekScalBlkMatSharedPtr
v_TransformedSchurCompl(int offset, const boost::shared_ptr<DNekScalBlkMat > &loc_mat);
virtual DNekScalMatSharedPtr
v_TransformedSchurCompl(int offset, const boost::shared_ptr<DNekScalMat > &loc_mat);
};
}
}
......
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