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

Conjugate gradient with low energy basis.

Get different R matrices from preconditioner.


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@4094 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent bc15701f
///////////////////////////////////////////////////////////////////////////////
//
// File GlobalLinSys.cpp
......@@ -309,6 +310,12 @@ namespace Nektar
return NullDNekMatSharedPtr;
}
void GlobalLinSys::v_InitObject()
{
NEKERROR(ErrorUtil::efatal,"Method does not exist" );
}
} //end of namespace
} //end of namespace
......@@ -87,6 +87,8 @@ namespace Nektar
const inline DNekMatSharedPtr &GetGmat(void) const;
inline void InitObject();
/// Solve the linear system for given input and output vectors
/// using a specified local to global map.
MULTI_REGIONS_EXPORT
......@@ -141,6 +143,8 @@ namespace Nektar
virtual const DNekMatSharedPtr& v_GetGmat(void) const;
virtual void v_InitObject();
static std::string lookupIds[];
static std::string def;
};
......@@ -197,6 +201,10 @@ namespace Nektar
return v_GetGmat();
}
inline void GlobalLinSys::InitObject()
{
v_InitObject();
}
} //end of namespace
} //end of namespace
......
......@@ -90,7 +90,7 @@ namespace Nektar
: GlobalLinSysIterative(pKey, pExpList, pLocToGloMap),
m_locToGloMap (pLocToGloMap)
{
ASSERTL1((pKey.GetGlobalSysSolnType()==eIterativeStaticCond)||
ASSERTL1((pKey.GetGlobalSysSolnType()==eIterativeStaticCond)||
(pKey.GetGlobalSysSolnType()==eIterativeMultiLevelStaticCond),
"This constructor is only valid when using static "
"condensation");
......@@ -100,10 +100,10 @@ namespace Nektar
"solution type");
// Allocate memory for top-level structure
SetupTopLevel(pLocToGloMap);
//SetupTopLevel(pLocToGloMap);
// Construct this level
Initialise(pLocToGloMap);
//Initialise(pLocToGloMap);
}
......@@ -131,6 +131,23 @@ namespace Nektar
}
void GlobalLinSysIterativeStaticCond::v_InitObject()
{
// Allocate memory for top-level structure
if (m_locToGloMap->GetPreconType() == MultiRegions::eLowEnergy)
{
SetupLowEnergyTopLevel(m_locToGloMap);
}
else
{
SetupTopLevel(m_locToGloMap);
}
// Construct this level
Initialise(m_locToGloMap);
}
/**
*
*/
......@@ -160,6 +177,7 @@ namespace Nektar
int nIntDofs = pLocToGloMap->GetNumGlobalCoeffs()
- nGlobBndDofs;
Array<OneD, NekDouble> F = m_wsp + nLocBndDofs;
if(nDirBndDofs && dirForcCalculated)
{
......@@ -172,6 +190,8 @@ namespace Nektar
NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,F+nDirBndDofs,
eWrapper);
NekVector<NekDouble> F_GlobBnd(nGlobBndDofs,F,eWrapper);
NekVector<NekDouble> F_LocBnd(nLocBndDofs,0.0);
NekVector<NekDouble> F_Int(nIntDofs,F+nGlobBndDofs,eWrapper);
NekVector<NekDouble> V_GlobBnd(nGlobBndDofs,out,eWrapper);
......@@ -184,39 +204,89 @@ namespace Nektar
if(nGlobHomBndDofs)
{
// 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;
if(pLocToGloMap->GetPreconType() != MultiRegions::eLowEnergy)
{
DNekScalBlkMat &BinvD = *m_BinvD;
DNekScalBlkMat &SchurCompl = *m_schurCompl;
// construct boundary forcing
if( nIntDofs && ((nDirBndDofs) && (!dirForcCalculated)
&& (atLastLevel)) )
{
//include dirichlet boundary forcing
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;
}
else
{
DNekScalBlkMat &BinvD = *m_BinvD;
V_LocBnd = BinvD*F_Int;
}
pLocToGloMap->AssembleBnd(V_LocBnd,V_GlobHomBndTmp,
nDirBndDofs);
F_HomBnd = F_HomBnd - V_GlobHomBndTmp;
DNekScalBlkMat &S1 = *m_S1Blk;
DNekScalBlkMat &R = *m_RBlk;
DNekScalBlkMat &BinvD = *m_BinvD;
if( nIntDofs && ((nDirBndDofs) && (!dirForcCalculated)
&& (atLastLevel)) )
{
pLocToGloMap->GlobalToLocalBnd(V_GlobBnd,V_LocBnd);
V_LocBnd = BinvD*F_Int+ S1*V_LocBnd;
}
else if((nDirBndDofs) && (!dirForcCalculated)
&& (atLastLevel))
{
pLocToGloMap->GlobalToLocalBnd(V_GlobBnd,V_LocBnd);
V_LocBnd = S1*V_LocBnd;
}
else
{
V_LocBnd = BinvD*F_Int;
}
pLocToGloMap->AssembleBnd(V_LocBnd,V_GlobHomBndTmp,
nDirBndDofs);
//F_bnd- B invD*F_int-S1*x
F_HomBnd = F_HomBnd - V_GlobHomBndTmp;
NekVector<NekDouble> fml(nLocBndDofs,0.0);
NekVector<NekDouble> fMultVector(nGlobBndDofs,1.0);
pLocToGloMap->GlobalToLocalBnd(fMultVector,fml);
pLocToGloMap->AssembleBnd(fml,fMultVector);
for(int i=0; i<nGlobBndDofs; ++i)
{
fMultVector[i]=1/fMultVector[i];
}
F_GlobBnd=F_GlobBnd*fMultVector;
pLocToGloMap->GlobalToLocalBnd(F_GlobBnd,F_LocBnd);
F_LocBnd=R*F_LocBnd;
pLocToGloMap->AssembleBnd(F_LocBnd,F_HomBnd, nDirBndDofs);
}
// solve boundary system
if(atLastLevel)
{
Array<OneD, NekDouble> offsetarray;
//Solve(F_HomBnd,V_GlobHomBnd);
//SolveLinearSystem(nGlobHomBndDofs, F+nDirBndDofs,offsetarray=out+nDirBndDofs);
Timer t;
t.Start();
......@@ -224,7 +294,31 @@ namespace Nektar
SolveLinearSystem(nGlobBndDofs, F, out, pLocToGloMap, nDirBndDofs);
t.Stop();
// std::cout << "time per solveLinearSystem = " << t.TimePerTest(1) << std::endl;
//transform back to original basis
if(pLocToGloMap->GetPreconType() == MultiRegions::eLowEnergy)
{
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;
}
}
else
{
......@@ -315,6 +409,7 @@ namespace Nektar
const Array<OneD,const unsigned int>& nint_size
= pLocToGloMap->GetNumLocalIntCoeffsPerPatch();
// Setup Block Matrix systems
MatrixStorage blkmatStorage = eDIAGONAL;
m_schurCompl = MemoryManager<DNekScalBlkMat>
......@@ -345,6 +440,80 @@ namespace Nektar
}
}
void GlobalLinSysIterativeStaticCond::SetupLowEnergyTopLevel(
const boost::shared_ptr<AssemblyMap>& pLocToGloMap)
{
int n;
int n_exp = m_expList.lock()->GetNumElmts();
const Array<OneD,const unsigned int>& nbdry_size
= pLocToGloMap->GetNumLocalBndCoeffsPerPatch();
const Array<OneD,const unsigned int>& nint_size
= pLocToGloMap->GetNumLocalIntCoeffsPerPatch();
v_UniqueMap();
m_precon = MemoryManager<Preconditioner>::AllocateSharedPtr(
GetSharedThisPtr(),m_locToGloMap);
// Setup Block Matrix systems
MatrixStorage blkmatStorage = eDIAGONAL;
m_schurCompl = MemoryManager<DNekScalBlkMat>
::AllocateSharedPtr(nbdry_size, nbdry_size, blkmatStorage);
m_BinvD = MemoryManager<DNekScalBlkMat>
::AllocateSharedPtr(nbdry_size, nint_size , blkmatStorage);
m_C = MemoryManager<DNekScalBlkMat>
::AllocateSharedPtr(nint_size , nbdry_size, blkmatStorage);
m_invD = MemoryManager<DNekScalBlkMat>
::AllocateSharedPtr(nint_size , nint_size , blkmatStorage);
//Variants of R matrices required for low energy preconditioning
m_RBlk = MemoryManager<DNekScalBlkMat>
::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
m_RTBlk = MemoryManager<DNekScalBlkMat>
::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
m_S1Blk = MemoryManager<DNekScalBlkMat>
::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
DNekMatSharedPtr m_R = m_precon->GetTransformationMatrix();
DNekMatSharedPtr m_RT = m_precon->GetTransposedTransformationMatrix();
for(n = 0; n < n_exp; ++n)
{
DNekScalBlkMatSharedPtr loc_mat = GetStaticCondBlock(m_expList.lock()->GetOffset_Elmt_Id(n));
DNekScalMatSharedPtr tmp_mat;
DNekScalMatSharedPtr m_S1=loc_mat->GetBlock(0,0);
DNekScalMat &S1 = (*m_S1);
int nRow=S1.GetRows();
NekDouble zero = 0.0;
NekDouble one = 1.0;
MatrixStorage storage = eFULL;
DNekMatSharedPtr m_S2 = MemoryManager<DNekMat>::AllocateSharedPtr(nRow,nRow,zero,storage);
DNekMatSharedPtr m_RS1 = MemoryManager<DNekMat>::AllocateSharedPtr(nRow,nRow,zero,storage);
//transformation matrices
DNekMat &R = (*m_R);
DNekMat &RT = (*m_RT);
//create low energy matrix
DNekMat &RS1 = (*m_RS1);
DNekMat &S2 = (*m_S2);
//setup S2
RS1=R*S1;
S2=RS1*RT;
m_schurCompl->SetBlock(n,n, tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,m_S2));
m_BinvD ->SetBlock(n,n, tmp_mat = loc_mat->GetBlock(0,1));
m_C ->SetBlock(n,n, tmp_mat = loc_mat->GetBlock(1,0));
m_invD ->SetBlock(n,n, tmp_mat = loc_mat->GetBlock(1,1));
m_S1Blk->SetBlock(n,n, tmp_mat = loc_mat->GetBlock(0,0));
m_RBlk->SetBlock(n,n, tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,m_R));
m_RTBlk->SetBlock(n,n, tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,m_RT));
}
}
/**
* Assemble the schur complement matrix from the block matrices stored
......@@ -676,7 +845,6 @@ namespace Nektar
else
{
// Do matrix multiply locally
NekVector<NekDouble> loc(nLocal, m_wsp, eWrapper);
m_locToGloMap->GlobalToLocalBnd(pInput, m_wsp);
loc = (*m_schurCompl)*loc;
......@@ -684,45 +852,10 @@ namespace Nektar
}
}
/**
* 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 GlobalLinSysIterativeStaticCond::v_ComputePreconditioner()
{
/*
ASSERTL1(m_gmat.get(),
"Matrix must be defined to compute preconditioner.");
ASSERTL1(!m_preconditioner.get(),
"Preconditioner has already been defined.");
int nGlobalBnd = m_locToGloMap->GetNumGlobalBndCoeffs();
int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
int n = m_gmat->GetRows();
m_map = m_locToGloMap->GetGlobalToUniversalBndMapUnique();
MatrixStorage storage = eDIAGONAL;
m_preconditioner = MemoryManager<DNekMat>::AllocateSharedPtr(n, n, storage);
DNekMat &M = (*m_preconditioner);
// Extract diagonal contributions
Array<OneD, NekDouble> vOutput(nGlobalBnd,0.0);
for (unsigned int i = 0; i < n; ++i)
{
vOutput[nDirBnd + i] = (*m_gmat)(i,i);
}
// Assemble diagonal contributions across processes
m_locToGloMap->UniversalAssembleBnd(vOutput);
// Populate preconditioner matrix
for (unsigned int i = 0; i < n; ++i)
{
M.SetValue(i,i,1.0/vOutput[nDirBnd + i]);
}
*/
}
//
}
void GlobalLinSysIterativeStaticCond::v_UniqueMap()
{
......
......@@ -63,8 +63,9 @@ namespace Nektar
const boost::shared_ptr<AssemblyMap>
&pLocToGloMap)
{
return MemoryManager<GlobalLinSysIterativeStaticCond>
::AllocateSharedPtr(pLinSysKey, pExpList, pLocToGloMap);
GlobalLinSysSharedPtr p = MemoryManager<GlobalLinSysIterativeStaticCond>::AllocateSharedPtr(pLinSysKey, pExpList, pLocToGloMap);
p->InitObject();
return p;
}
/// Name of class
......@@ -101,6 +102,11 @@ namespace Nektar
DNekScalBlkMatSharedPtr m_C;
DNekScalBlkMatSharedPtr m_invD;
// Block matrices for low energy
DNekScalBlkMatSharedPtr m_RBlk;
DNekScalBlkMatSharedPtr m_RTBlk;
DNekScalBlkMatSharedPtr m_S1Blk;
/// Globally assembled Schur complement matrix at this level
GlobalMatrixSharedPtr m_globalSchurCompl;
......@@ -110,6 +116,8 @@ namespace Nektar
// Workspace array for matrix multiplication
Array<OneD, NekDouble> m_wsp;
PreconditionerSharedPtr m_precon;
/// Solve the linear system for given input and output vectors
/// using a specified local to global map.
virtual void v_Solve(
......@@ -119,6 +127,8 @@ namespace Nektar
const Array<OneD, const NekDouble> &dirForcing
= NullNekDouble1DArray);
virtual void v_InitObject();
/// Initialise this object
void Initialise(
const boost::shared_ptr<AssemblyMap>& locToGloMap);
......@@ -128,6 +138,9 @@ namespace Nektar
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);
......
......@@ -114,7 +114,8 @@ namespace Nektar
eInverseLinear,
eLowEnergy,
eLinearLowEnergy,
eBlock
eBlock,
eLocalLowEnergy
};
const char* const PreconditionerTypeMap[] =
......@@ -124,7 +125,8 @@ namespace Nektar
"InverseLinear",
"LowEnergy",
"LinearLowEnergy",
"Block"
"Block",
"LocalLowEnergy"
};
......
This diff is collapsed.
......@@ -62,6 +62,8 @@ namespace Nektar
const DNekMatSharedPtr& GetTransposedTransformationMatrix() const;
const DNekMatSharedPtr& GetInverseTransformationMatrix() const;
const DNekMatSharedPtr& GetInverseTransposedTransformationMatrix() const;
protected:
......@@ -81,6 +83,7 @@ namespace Nektar
DNekMatSharedPtr m_transformationmatrix;
DNekMatSharedPtr m_inversetransformationmatrix;
DNekMatSharedPtr m_transposedtransformationmatrix;
DNekMatSharedPtr m_inversetransposedtransformationmatrix;
DNekMatSharedPtr m_efedgefacecoupling;
DNekMatSharedPtr m_effacefacecoupling;
DNekMatSharedPtr m_edgefacetransformmatrix;
......
Supports Markdown
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