Commit 8869223e authored by David Moxey's avatar David Moxey
Browse files

Working iterative static cond without low energy precon

parent c0e7ecfb
......@@ -129,131 +129,6 @@ namespace Nektar
}
#if 0
/**
*
*/
void GlobalLinSysDirectStaticCond::v_Solve(
const Array<OneD, const NekDouble> &in,
Array<OneD, NekDouble> &out,
const AssemblyMapSharedPtr &pLocToGloMap,
const Array<OneD, const NekDouble> &dirForcing)
{
bool dirForcCalculated = (bool) dirForcing.num_elements();
bool atLastLevel = pLocToGloMap->AtLastLevel();
int nGlobDofs = pLocToGloMap->GetNumGlobalCoeffs();
int nGlobBndDofs = pLocToGloMap->GetNumGlobalBndCoeffs();
int nDirBndDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
int nLocBndDofs = pLocToGloMap->GetNumLocalBndCoeffs();
int nIntDofs = pLocToGloMap->GetNumGlobalCoeffs()
- nGlobBndDofs;
Array<OneD, NekDouble> F(nGlobDofs);
Array<OneD, NekDouble> tmp;
if(nDirBndDofs && dirForcCalculated)
{
Vmath::Vsub(nGlobDofs,in.get(),1,dirForcing.get(),1,F.get(),1);
}
else
{
Vmath::Vcopy(nGlobDofs,in.get(),1,F.get(),1);
}
NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,tmp=F+nDirBndDofs,
eWrapper);
NekVector<NekDouble> F_Int(nIntDofs,tmp=F+nGlobBndDofs,eWrapper);
NekVector<NekDouble> V_GlobBnd(nGlobBndDofs,out,eWrapper);
NekVector<NekDouble> V_GlobHomBnd(nGlobHomBndDofs,
tmp=out+nDirBndDofs,
eWrapper);
NekVector<NekDouble> V_Int(nIntDofs,tmp=out+nGlobBndDofs,eWrapper);
NekVector<NekDouble> V_LocBnd(nLocBndDofs,0.0);
NekVector<NekDouble> V_GlobHomBndTmp(nGlobHomBndDofs,0.0);
//zero GlobHomBnd so that we ensure we are solving for
//full problem rather than perturbation from initial
//condition in this case
Vmath::Zero(nGlobHomBndDofs,tmp = out+nDirBndDofs,1);
if(nGlobHomBndDofs)
{
if(nIntDofs || ((nDirBndDofs) && (!dirForcCalculated)
&& (atLastLevel)) )
{
// construct boundary forcing
if( nIntDofs && ((nDirBndDofs) && (!dirForcCalculated)
&& (atLastLevel)) )
{
//include dirichlet boundary forcing
DNekScalBlkMat &BinvD = *m_BinvD;
DNekScalBlkMat &SchurCompl = *m_schurCompl;
pLocToGloMap->GlobalToLocalBnd(V_GlobBnd,V_LocBnd);
V_LocBnd = BinvD*F_Int + SchurCompl*V_LocBnd;
}
else if((nDirBndDofs) && (!dirForcCalculated)
&& (atLastLevel))
{
//include dirichlet boundary forcing
DNekScalBlkMat &SchurCompl = *m_schurCompl;
pLocToGloMap->GlobalToLocalBnd(V_GlobBnd,V_LocBnd);
V_LocBnd = SchurCompl*V_LocBnd;
}
else
{
DNekScalBlkMat &BinvD = *m_BinvD;
V_LocBnd = BinvD*F_Int;
}
pLocToGloMap->AssembleBnd(V_LocBnd,V_GlobHomBndTmp,
nDirBndDofs);
F_HomBnd = F_HomBnd - V_GlobHomBndTmp;
}
// solve boundary system
if(atLastLevel)
{
m_linSys->Solve(F_HomBnd,V_GlobHomBnd);
}
else
{
m_recursiveSchurCompl->Solve(F,
V_GlobBnd.GetPtr(),
pLocToGloMap->GetNextLevelLocalToGlobalMap());
}
}
// solve interior system
if(nIntDofs)
{
DNekScalBlkMat &invD = *m_invD;
if(nGlobHomBndDofs || nDirBndDofs)
{
DNekScalBlkMat &C = *m_C;
if(dirForcCalculated && nDirBndDofs)
{
pLocToGloMap->GlobalToLocalBnd(V_GlobHomBnd,V_LocBnd,
nDirBndDofs);
}
else
{
pLocToGloMap->GlobalToLocalBnd(V_GlobBnd,V_LocBnd);
}
F_Int = F_Int - C*V_LocBnd;
}
V_Int = invD*F_Int;
}
}
#endif
MatrixStorage GlobalLinSysDirectStaticCond::DetermineMatrixStorage(
const AssemblyMapSharedPtr &pLocToGloMap)
......
......@@ -95,27 +95,9 @@ namespace Nektar
const AssemblyMapSharedPtr &plocToGloMap,
const int nDir)
{
// Get vector sizes
int nNonDir = nGlobal - nDir;
Array<OneD, NekDouble> tmp;
if (m_useProjection)
{
DoAconjugateProjection(nGlobal, pInput, pOutput, plocToGloMap, nDir);
if (0)
{
// check correctness: solve the same system with plain CG and compare
Array<OneD, NekDouble> cg_s (nGlobal, 0.0);
NekVector<NekDouble> cg (nNonDir, tmp = cg_s + nDir, eWrapper);
NekVector<NekDouble> x (nNonDir, tmp = pOutput + nDir, eWrapper);
DoConjugateGradient(nGlobal, pInput, cg_s, plocToGloMap, nDir);
cg -= x;
NekDouble norm = CalculateAnorm(nGlobal, cg_s, nDir);
std::cout << "norm of solutions difference is = " << norm << std::endl;
}
}
else
{
......@@ -130,7 +112,6 @@ namespace Nektar
* in order to speed up successive linear solves with
* right-hand sides arising from time-dependent discretisations.
* (P.F.Fischer, Comput. Methods Appl. Mech. Engrg. 163, 1998)
*
*/
void GlobalLinSysIterative::DoAconjugateProjection(
const int nGlobal,
......@@ -372,19 +353,22 @@ namespace Nektar
* @param pOutput Solution vector of all DOFs.  
*/
void GlobalLinSysIterative::DoConjugateGradient(
const int nGlobal,
const Array<OneD,const NekDouble> &pInput,
Array<OneD, NekDouble> &pOutput,
const AssemblyMapSharedPtr &plocToGloMap,
const int nDir)
const int nGlobal,
const Array<OneD,const NekDouble> &pInput,
Array<OneD, NekDouble> &pOutput,
const AssemblyMapSharedPtr &plocToGloMap,
const int nDir)
{
if (!m_precon)
{
MultiRegions::PreconditionerType pType = plocToGloMap->GetPreconType();
std::string PreconType = MultiRegions::PreconditionerTypeMap[pType];
MultiRegions::PreconditionerType pType
= plocToGloMap->GetPreconType();
std::string PreconType
= MultiRegions::PreconditionerTypeMap[pType];
v_UniqueMap();
m_precon = GetPreconFactory().CreateInstance(PreconType,GetSharedThisPtr(),plocToGloMap);
m_precon -> BuildPreconditioner();
m_precon = GetPreconFactory().CreateInstance(
PreconType,GetSharedThisPtr(),plocToGloMap);
m_precon->BuildPreconditioner();
}
// Get the communicator for performing data exchanges
......
......@@ -49,15 +49,14 @@ namespace Nektar
class ExpList;
/// A global linear system.
class GlobalLinSysIterative : public GlobalLinSys
class GlobalLinSysIterative : virtual public GlobalLinSys
{
public:
/// Constructor for full direct matrix solve.
MULTI_REGIONS_EXPORT GlobalLinSysIterative(
const GlobalLinSysKey &pKey,
const boost::weak_ptr<ExpList> &pExpList,
const boost::shared_ptr<AssemblyMap>
&pLocToGloMap);
const GlobalLinSysKey &pKey,
const boost::weak_ptr<ExpList> &pExpList,
const boost::shared_ptr<AssemblyMap> &pLocToGloMap);
MULTI_REGIONS_EXPORT virtual ~GlobalLinSysIterative();
......@@ -74,7 +73,6 @@ namespace Nektar
PreconditionerSharedPtr m_precon;
MultiRegions::PreconditionerType m_precontype;
PreconditionerSharedPtr m_lowEnergyPrecon;
int m_totalIterations;
......@@ -91,7 +89,6 @@ namespace Nektar
/// Total counter of previous solutions
int m_numPrevSols;
/// A-conjugate projection technique
void DoAconjugateProjection(
const int pNumRows,
......@@ -110,16 +107,10 @@ namespace Nektar
void Set_Rhs_Magnitude(const NekVector<NekDouble> &pIn);
virtual void v_UniqueMap() = 0;
private:
void printArray(
const std::string& msg,
const Array<OneD, const NekDouble> &in,
const int len,
const int offset);
void UpdateKnownSolutions(
const int pGlobalBndDofs,
const Array<OneD,const NekDouble> &pSolution,
......@@ -142,8 +133,6 @@ namespace Nektar
virtual void v_DoMatrixMultiply(
const Array<OneD, NekDouble>& pInput,
Array<OneD, NekDouble>& pOutput) = 0;
virtual void v_UniqueMap() = 0;
};
}
}
......
......@@ -70,7 +70,8 @@ namespace Nektar
const GlobalLinSysKey &pKey,
const boost::weak_ptr<ExpList> &pExp,
const boost::shared_ptr<AssemblyMap> &pLocToGloMap)
: GlobalLinSysIterative(pKey, pExp, pLocToGloMap)
: GlobalLinSys (pKey, pExp, pLocToGloMap),
GlobalLinSysIterative(pKey, pExp, pLocToGloMap)
{
ASSERTL1(m_linSysKey.GetGlobalSysSolnType()==eIterativeFull,
"This routine should only be used when using an Iterative "
......
......@@ -37,6 +37,7 @@
#include <MultiRegions/GlobalMatrix.h>
#include <MultiRegions/GlobalLinSysIterative.h>
#include <MultiRegions/GlobalLinSysStaticCond.h>
#include <LibUtilities/LinearAlgebra/SparseMatrixFwd.hpp>
......@@ -68,7 +69,8 @@ namespace Nektar
/// A global linear system.
class GlobalLinSysIterativeStaticCond : public GlobalLinSysIterative
class GlobalLinSysIterativeStaticCond : virtual public GlobalLinSysIterative,
virtual public GlobalLinSysStaticCond
{
public:
typedef NekSparseDiagBlkMatrix<StorageSmvBsr<NekDouble> >
......@@ -113,25 +115,18 @@ namespace Nektar
virtual ~GlobalLinSysIterativeStaticCond();
protected:
virtual int v_GetNumBlocks();
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,
const DNekScalBlkMatSharedPtr pSchurCompl,
const DNekScalBlkMatSharedPtr pBinvD,
const DNekScalBlkMatSharedPtr pC,
const DNekScalBlkMatSharedPtr pInvD,
const boost::shared_ptr<AssemblyMap> &locToGloMap);
private:
/// Schur complement for Direct Static Condensation.
GlobalLinSysIterativeStaticCondSharedPtr m_recursiveSchurCompl;
/// Block Schur complement matrix.
DNekScalBlkMatSharedPtr m_schurCompl;
/// Block \f$ BD^{-1} \f$ matrix.
DNekScalBlkMatSharedPtr m_BinvD;
/// Block \f$ C \f$ matrix.
DNekScalBlkMatSharedPtr m_C;
/// Block \f$ D^{-1} \f$ matrix.
DNekScalBlkMatSharedPtr m_invD;
/// Block matrices for low energy
DNekScalBlkMatSharedPtr m_RBlk;
DNekScalBlkMatSharedPtr m_RTBlk;
DNekScalBlkMatSharedPtr m_S1Blk;
/// Dense storage for block Schur complement matrix
std::vector<double> m_storage;
/// Vector of pointers to local matrix data
......@@ -140,58 +135,22 @@ namespace Nektar
Array<OneD, unsigned int> m_rows;
/// Scaling factors for local matrices
Array<OneD, NekDouble> m_scale;
/// Sparse representation of Schur complement matrix at this level
DNekSmvBsrDiagBlkMatSharedPtr m_sparseSchurCompl;
/// Local to global map.
boost::shared_ptr<AssemblyMap> m_locToGloMap;
/// Workspace array for matrix multiplication
Array<OneD, NekDouble> m_wsp;
/// Preconditioner object.
PreconditionerSharedPtr m_precon;
/// Utility strings
static std::string storagedef;
static std::string storagelookupIds[];
/// Solve the linear system for given input and output vectors
/// using a specified local to global map.
virtual void v_Solve(
const Array<OneD, const NekDouble> &in,
Array<OneD, NekDouble> &out,
const AssemblyMapSharedPtr &locToGloMap,
const Array<OneD, const NekDouble> &dirForcing
= NullNekDouble1DArray);
virtual void v_InitObject();
/// Initialise this object
void Initialise(
const boost::shared_ptr<AssemblyMap>& locToGloMap);
/// Set up the storage for the Schur complement or the top level
/// of the multi-level Schur complement.
void SetupTopLevel(
const boost::shared_ptr<AssemblyMap>& locToGloMap);
void SetupLowEnergyTopLevel(
const boost::shared_ptr<AssemblyMap>& locToGloMap);
/// Assemble the Schur complement matrix.
void AssembleSchurComplement(
const boost::shared_ptr<AssemblyMap>& locToGloMap);
void v_AssembleSchurComplement(
const boost::shared_ptr<AssemblyMap> locToGloMap);
/// Prepares local representation of Schur complement
/// stored as a sparse block-diagonal matrix.
void PrepareLocalSchurComplement();
///
void ConstructNextLevelCondensedSystem(
const boost::shared_ptr<AssemblyMap>& locToGloMap);
/// Perform a Shur-complement matrix multiply operation.
virtual void v_DoMatrixMultiply(
const Array<OneD, NekDouble>& pInput,
......
......@@ -92,10 +92,6 @@ namespace Nektar
}
virtual int v_GetNumBlocks();
virtual DNekScalBlkMatSharedPtr v_GetStaticCondBlock(unsigned int n)
{
return GlobalLinSys::GetStaticCondBlock(n);
}
virtual GlobalLinSysStaticCondSharedPtr v_Recurse(
const GlobalLinSysKey &mkey,
......
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