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

Tidy up linear space preconditioner option to use either PETSc or Xxt

parent f882c192
...@@ -11,6 +11,7 @@ IF( NEKTAR_USE_PETSC ) ...@@ -11,6 +11,7 @@ IF( NEKTAR_USE_PETSC )
ELSE (THIRDPARTY_BUILD_PETSC) ELSE (THIRDPARTY_BUILD_PETSC)
INCLUDE (FindPETSc) INCLUDE (FindPETSc)
MESSAGE(STATUS "Found PETSc: ${PETSC_LIBRARIES}") MESSAGE(STATUS "Found PETSc: ${PETSC_LIBRARIES}")
ADD_DEFINITIONS(-DNEKTAR_USING_PETSC)
ENDIF (THIRDPARTY_BUILD_PETSC) ENDIF (THIRDPARTY_BUILD_PETSC)
ENDIF( NEKTAR_USE_PETSC ) ENDIF( NEKTAR_USE_PETSC )
......
...@@ -309,7 +309,8 @@ namespace Nektar ...@@ -309,7 +309,8 @@ namespace Nektar
m_solnType = solnTypeOld; m_solnType = solnTypeOld;
ASSERTL1(m_solnType==eDirectMultiLevelStaticCond ASSERTL1(m_solnType==eDirectMultiLevelStaticCond
||m_solnType==eIterativeMultiLevelStaticCond ||m_solnType==eIterativeMultiLevelStaticCond
||m_solnType==eXxtMultiLevelStaticCond, ||m_solnType==eXxtMultiLevelStaticCond
||m_solnType==ePETScMultiLevelStaticCond,
"This method should only be called for in " "This method should only be called for in "
"case of multi-level static condensation."); "case of multi-level static condensation.");
m_staticCondLevel = newLevel; m_staticCondLevel = newLevel;
...@@ -639,7 +640,8 @@ namespace Nektar ...@@ -639,7 +640,8 @@ namespace Nektar
return result; return result;
} }
boost::shared_ptr<AssemblyMap> AssemblyMap::v_XxtLinearSpaceMap(const ExpList &locexp) boost::shared_ptr<AssemblyMap> AssemblyMap::v_LinearSpaceMap(
const ExpList &locexp, GlobalSysSolnType solnType)
{ {
ASSERTL0(false, "Not defined for this sub class"); ASSERTL0(false, "Not defined for this sub class");
static boost::shared_ptr<AssemblyMap> result; static boost::shared_ptr<AssemblyMap> result;
...@@ -802,9 +804,9 @@ namespace Nektar ...@@ -802,9 +804,9 @@ namespace Nektar
return v_GetExtraDirEdges(); return v_GetExtraDirEdges();
} }
boost::shared_ptr<AssemblyMap> AssemblyMap::XxtLinearSpaceMap(const ExpList &locexp) boost::shared_ptr<AssemblyMap> AssemblyMap::LinearSpaceMap(const ExpList &locexp, GlobalSysSolnType solnType)
{ {
return v_XxtLinearSpaceMap(locexp); return v_LinearSpaceMap(locexp, solnType);
} }
int AssemblyMap::GetLocalToGlobalBndMap(const int i) const int AssemblyMap::GetLocalToGlobalBndMap(const int i) const
...@@ -1210,7 +1212,6 @@ namespace Nektar ...@@ -1210,7 +1212,6 @@ namespace Nektar
= m_session->GetComm()->GetRowComm(); = m_session->GetComm()->GetRowComm();
bool isRoot = vRowComm->GetRank() == 0; bool isRoot = vRowComm->GetRank() == 0;
int n = vRowComm->GetSize(); int n = vRowComm->GetSize();
int p = vRowComm->GetRank();
int i; int i;
// Determine number of global degrees of freedom. // Determine number of global degrees of freedom.
......
...@@ -257,7 +257,7 @@ namespace Nektar ...@@ -257,7 +257,7 @@ namespace Nektar
MULTI_REGIONS_EXPORT const Array<OneD, const int>& MULTI_REGIONS_EXPORT const Array<OneD, const int>&
GetExtraDirEdges(); GetExtraDirEdges();
MULTI_REGIONS_EXPORT boost::shared_ptr<AssemblyMap> XxtLinearSpaceMap(const ExpList &locexp); MULTI_REGIONS_EXPORT boost::shared_ptr<AssemblyMap> LinearSpaceMap(const ExpList &locexp, GlobalSysSolnType solnType);
/// Returns the bandwidth of the boundary system. /// Returns the bandwidth of the boundary system.
MULTI_REGIONS_EXPORT int GetBndSystemBandWidth() const; MULTI_REGIONS_EXPORT int GetBndSystemBandWidth() const;
...@@ -473,8 +473,8 @@ namespace Nektar ...@@ -473,8 +473,8 @@ namespace Nektar
v_GetExtraDirEdges(); v_GetExtraDirEdges();
/// Generate a linear space mapping from existing mapping /// Generate a linear space mapping from existing mapping
virtual boost::shared_ptr<AssemblyMap> v_XxtLinearSpaceMap virtual boost::shared_ptr<AssemblyMap> v_LinearSpaceMap(
(const ExpList &locexp); const ExpList &locexp, GlobalSysSolnType solnType);
}; };
......
...@@ -351,11 +351,12 @@ namespace Nektar ...@@ -351,11 +351,12 @@ namespace Nektar
* @brief Construct an AssemblyMapCG object which corresponds to the * @brief Construct an AssemblyMapCG object which corresponds to the
* linear space of the current object. * linear space of the current object.
* *
* This function is used in an XXT solve to apply a linear space * This function is used to create a linear-space assembly map, which is
* preconditioner in the conjugate gradient solve. * then used in the linear space preconditioner in the conjugate
* gradient solve.
*/ */
AssemblyMapSharedPtr AssemblyMapCG::v_XxtLinearSpaceMap( AssemblyMapSharedPtr AssemblyMapCG::v_LinearSpaceMap(
const ExpList &locexp) const ExpList &locexp, GlobalSysSolnType solnType)
{ {
AssemblyMapCGSharedPtr returnval; AssemblyMapCGSharedPtr returnval;
...@@ -368,7 +369,7 @@ namespace Nektar ...@@ -368,7 +369,7 @@ namespace Nektar
// Get Default Map and turn off any searched values. // Get Default Map and turn off any searched values.
returnval = MemoryManager<AssemblyMapCG> returnval = MemoryManager<AssemblyMapCG>
::AllocateSharedPtr(m_session); ::AllocateSharedPtr(m_session);
returnval->m_solnType = ePETScFullMatrix; //XxtFullMatrix; returnval->m_solnType = solnType;
returnval->m_preconType = eNull; returnval->m_preconType = eNull;
returnval->m_maxStaticCondLevel = 0; returnval->m_maxStaticCondLevel = 0;
returnval->m_signChange = false; returnval->m_signChange = false;
......
...@@ -195,7 +195,8 @@ namespace Nektar ...@@ -195,7 +195,8 @@ namespace Nektar
MULTI_REGIONS_EXPORT virtual const Array<OneD, const int>& v_GetExtraDirEdges(); MULTI_REGIONS_EXPORT virtual const Array<OneD, const int>& v_GetExtraDirEdges();
MULTI_REGIONS_EXPORT virtual AssemblyMapSharedPtr v_XxtLinearSpaceMap(const ExpList &locexp); MULTI_REGIONS_EXPORT virtual AssemblyMapSharedPtr v_LinearSpaceMap(
const ExpList &locexp, GlobalSysSolnType solnType);
}; };
......
...@@ -123,13 +123,11 @@ IF(NEKTAR_USE_PETSC) ...@@ -123,13 +123,11 @@ IF(NEKTAR_USE_PETSC)
GlobalLinSysPETSc.h GlobalLinSysPETSc.h
GlobalLinSysPETScFull.h GlobalLinSysPETScFull.h
GlobalLinSysPETScStaticCond.h GlobalLinSysPETScStaticCond.h
PreconditionerPETScLinear.h
) )
SET(MULTI_REGIONS_SOURCES ${MULTI_REGIONS_SOURCES} SET(MULTI_REGIONS_SOURCES ${MULTI_REGIONS_SOURCES}
GlobalLinSysPETSc.cpp GlobalLinSysPETSc.cpp
GlobalLinSysPETScFull.cpp GlobalLinSysPETScFull.cpp
GlobalLinSysPETScStaticCond.cpp GlobalLinSysPETScStaticCond.cpp
PreconditionerPETScLinear.cpp
) )
ENDIF(NEKTAR_USE_PETSC) ENDIF(NEKTAR_USE_PETSC)
......
...@@ -127,8 +127,7 @@ namespace Nektar ...@@ -127,8 +127,7 @@ namespace Nektar
eLowEnergy, eLowEnergy,
eLinearWithLowEnergy, eLinearWithLowEnergy,
eBlock, eBlock,
eLinearWithBlock, eLinearWithBlock
eLinearPETSc
}; };
const char* const PreconditionerTypeMap[] = const char* const PreconditionerTypeMap[] =
...@@ -140,8 +139,7 @@ namespace Nektar ...@@ -140,8 +139,7 @@ namespace Nektar
"LowEnergyBlock", "LowEnergyBlock",
"FullLinearSpaceWithLowEnergyBlock", "FullLinearSpaceWithLowEnergyBlock",
"Block", "Block",
"FullLinearSpaceWithBlock", "FullLinearSpaceWithBlock"
"LinearPETSc"
}; };
......
...@@ -42,7 +42,7 @@ namespace Nektar ...@@ -42,7 +42,7 @@ namespace Nektar
{ {
namespace MultiRegions namespace MultiRegions
{ {
std::string Preconditioner::lookupIds[9] = { std::string Preconditioner::lookupIds[8] = {
LibUtilities::SessionReader::RegisterEnumValue( LibUtilities::SessionReader::RegisterEnumValue(
"Preconditioner", "Null", eNull), "Preconditioner", "Null", eNull),
LibUtilities::SessionReader::RegisterEnumValue( LibUtilities::SessionReader::RegisterEnumValue(
...@@ -61,8 +61,6 @@ namespace Nektar ...@@ -61,8 +61,6 @@ namespace Nektar
"Preconditioner", "Block",eBlock), "Preconditioner", "Block",eBlock),
LibUtilities::SessionReader::RegisterEnumValue( LibUtilities::SessionReader::RegisterEnumValue(
"Preconditioner", "FullLinearSpaceWithBlock",eLinearWithBlock), "Preconditioner", "FullLinearSpaceWithBlock",eLinearWithBlock),
LibUtilities::SessionReader::RegisterEnumValue(
"Preconditioner", "LinearPETSc",eLinearPETSc),
}; };
std::string Preconditioner::def = std::string Preconditioner::def =
LibUtilities::SessionReader::RegisterDefaultSolverInfo( LibUtilities::SessionReader::RegisterDefaultSolverInfo(
......
...@@ -39,6 +39,11 @@ ...@@ -39,6 +39,11 @@
#include <MultiRegions/GlobalLinSysIterativeStaticCond.h> #include <MultiRegions/GlobalLinSysIterativeStaticCond.h>
#include <MultiRegions/GlobalLinSys.h> #include <MultiRegions/GlobalLinSys.h>
#include <MultiRegions/GlobalLinSysXxtFull.h> #include <MultiRegions/GlobalLinSysXxtFull.h>
#ifdef NEKTAR_USING_PETSC
#include <MultiRegions/GlobalLinSysPETScFull.h>
#endif
#include <LocalRegions/MatrixKey.h> #include <LocalRegions/MatrixKey.h>
#include <math.h> #include <math.h>
...@@ -56,7 +61,22 @@ namespace Nektar ...@@ -56,7 +61,22 @@ namespace Nektar
PreconditionerLinear::create, PreconditionerLinear::create,
"Full Linear space inverse Preconditioning"); "Full Linear space inverse Preconditioning");
/** std::string PreconditionerLinear::solveType =
LibUtilities::SessionReader::RegisterDefaultSolverInfo(
"LinearPreconSolver",
"Xxt");
std::string PreconditionerLinear::solveTypeIds[] = {
LibUtilities::SessionReader::RegisterEnumValue(
"LinearPreconSolver",
"PETSc",
MultiRegions::eLinearPreconPETSc),
LibUtilities::SessionReader::RegisterEnumValue(
"LinearPreconSolver",
"Xxt",
MultiRegions::eLinearPreconXxt)
};
/**
* @class PreconditionerLinear * @class PreconditionerLinear
* *
* This class implements preconditioning for the conjugate * This class implements preconditioning for the conjugate
...@@ -76,36 +96,67 @@ namespace Nektar ...@@ -76,36 +96,67 @@ namespace Nektar
void PreconditionerLinear::v_BuildPreconditioner() void PreconditionerLinear::v_BuildPreconditioner()
{ {
GlobalSysSolnType solvertype=m_locToGloMap->GetGlobalSysSolnType(); GlobalSysSolnType sType = m_locToGloMap->GetGlobalSysSolnType();
if(solvertype != eIterativeStaticCond) ASSERTL0(sType == eIterativeStaticCond,
{ "This type of preconditioning is not implemented "
ASSERTL0(0,"This type of preconditioning is not implemented for this solver"); "for this solver");
}
boost::shared_ptr<MultiRegions::ExpList> boost::shared_ptr<MultiRegions::ExpList>
expList=((m_linsys.lock())->GetLocMat()).lock(); expList=((m_linsys.lock())->GetLocMat()).lock();
m_vertLocToGloMap = m_locToGloMap->XxtLinearSpaceMap(*expList);
// Generate XXT system. LinearPreconSolver solveType =
if(m_linsys.lock()->GetKey().GetMatrixType() == StdRegions::eMass) expList->GetSession()->GetSolverInfoAsEnum<LinearPreconSolver>(
{ "LinearPreconSolver");
GlobalLinSysKey preconKey(StdRegions::ePreconLinearSpaceMass,
m_vertLocToGloMap); GlobalSysSolnType linSolveType;
m_vertLinsys = MemoryManager<GlobalLinSysXxtFull>::
AllocateSharedPtr(preconKey,expList,m_vertLocToGloMap); switch(solveType)
}
else
{ {
GlobalLinSysKey preconKey(StdRegions::ePreconLinearSpace, case eLinearPreconXxt:
m_vertLocToGloMap, {
(m_linsys.lock())->GetKey().GetConstFactors()); linSolveType = eXxtFullMatrix;
m_vertLinsys = MemoryManager<GlobalLinSysXxtFull>:: break;
AllocateSharedPtr(preconKey,expList,m_vertLocToGloMap); }
case eLinearPreconPETSc:
{
#ifdef NEKTAR_USING_PETSC
linSolveType = ePETScFullMatrix;
#else
ASSERTL0(false, "Nektar++ has not been compiled with "
"PETSc support.");
#endif
}
} }
m_vertLocToGloMap = m_locToGloMap->LinearSpaceMap(*expList, linSolveType);
// Generate linear solve system.
StdRegions::MatrixType mType =
m_linsys.lock()->GetKey().GetMatrixType() == StdRegions::eMass ?
StdRegions::ePreconLinearSpaceMass :
StdRegions::ePreconLinearSpace;
GlobalLinSysKey preconKey(mType, m_vertLocToGloMap);
switch(solveType)
{
case eLinearPreconXxt:
{
m_vertLinsys = MemoryManager<GlobalLinSysXxtFull>::
AllocateSharedPtr(preconKey,expList,m_vertLocToGloMap);
break;
}
case eLinearPreconPETSc:
{
#ifdef NEKTAR_USING_PETSC
m_vertLinsys = MemoryManager<GlobalLinSysPETScFull>::
AllocateSharedPtr(preconKey,expList,m_vertLocToGloMap);
#else
ASSERTL0(false, "Nektar++ has not been compiled with "
"PETSc support.");
#endif
}
}
} }
/** /**
......
...@@ -46,6 +46,12 @@ namespace Nektar ...@@ -46,6 +46,12 @@ namespace Nektar
{ {
namespace MultiRegions namespace MultiRegions
{ {
enum LinearPreconSolver
{
eLinearPreconXxt,
eLinearPreconPETSc
};
class PreconditionerLinear; class PreconditionerLinear;
typedef boost::shared_ptr<PreconditionerLinear> PreconditionerLinearSharedPtr; typedef boost::shared_ptr<PreconditionerLinear> PreconditionerLinearSharedPtr;
...@@ -78,6 +84,8 @@ namespace Nektar ...@@ -78,6 +84,8 @@ namespace Nektar
boost::shared_ptr<AssemblyMap> m_vertLocToGloMap; boost::shared_ptr<AssemblyMap> m_vertLocToGloMap;
private: private:
static std::string solveType;
static std::string solveTypeIds[];
virtual void v_InitObject(); virtual void v_InitObject();
......
///////////////////////////////////////////////////////////////////////////////
//
// File PreconditionerPETScLinear.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: PreconditionerPETScLinear definition
//
///////////////////////////////////////////////////////////////////////////////
#include <LibUtilities/BasicUtils/VDmathArray.hpp>
#include <MultiRegions/PreconditionerPETScLinear.h>
#include <MultiRegions/GlobalMatrixKey.h>
#include <MultiRegions/GlobalLinSysIterativeStaticCond.h>
#include <MultiRegions/GlobalLinSys.h>
#include <MultiRegions/GlobalLinSysPETScFull.h>
#include <LocalRegions/MatrixKey.h>
#include <math.h>
namespace Nektar
{
namespace MultiRegions
{
/**
* Registers the class with the Factory.
*/
string PreconditionerPETScLinear::className1
= GetPreconFactory().RegisterCreatorFunction(
"LinearPETSc",
PreconditionerPETScLinear::create,
"PETSc Preconditioning");
/**
* @class PreconditionerPETScLinear
*
* This class implements preconditioning for the conjugate
* gradient matrix solver.
*/
PreconditionerPETScLinear::PreconditionerPETScLinear(
const boost::shared_ptr<GlobalLinSys> &plinsys,
const AssemblyMapSharedPtr &pLocToGloMap)
: Preconditioner(plinsys, pLocToGloMap)
{
}
void PreconditionerPETScLinear::v_InitObject()
{
}
void PreconditionerPETScLinear::v_BuildPreconditioner()
{
GlobalSysSolnType solvertype=m_locToGloMap->GetGlobalSysSolnType();
boost::shared_ptr<MultiRegions::ExpList>
expList=((m_linsys.lock())->GetLocMat()).lock();
m_vertLocToGloMap = m_locToGloMap->XxtLinearSpaceMap(*expList);
// Generate PETSc system.
GlobalLinSysKey preconKey(StdRegions::ePreconLinearSpace,
m_vertLocToGloMap,
(m_linsys.lock())->GetKey().GetConstFactors());
m_vertLinsys = MemoryManager<GlobalLinSysPETScFull>::
AllocateSharedPtr(preconKey,expList,m_vertLocToGloMap);
}
/**
*
*/
void PreconditionerPETScLinear::v_DoPreconditioner(
const Array<OneD, NekDouble>& pInput,
Array<OneD, NekDouble>& pOutput)
{
v_DoPreconditionerWithNonVertOutput(pInput,pOutput,NullNekDouble1DArray);
}
/**
*
*/
void PreconditionerPETScLinear::v_DoPreconditionerWithNonVertOutput(
const Array<OneD, NekDouble>& pInput,
Array<OneD, NekDouble>& pOutput,
const Array<OneD, NekDouble>& pNonVertOutput)
{
GlobalSysSolnType solvertype=m_locToGloMap->GetGlobalSysSolnType();
switch(solvertype)
{
case MultiRegions::eIterativeStaticCond:
{
int i,val;
int nloc = m_vertLocToGloMap->GetNumLocalCoeffs();
int nglo = m_vertLocToGloMap->GetNumGlobalCoeffs();
// mapping from full space to vertices
Array<OneD, int> LocToGloBnd = m_vertLocToGloMap->GetLocalToGlobalBndMap();
// Global to local for linear solver (different from above)
Array<OneD, int> LocToGlo = m_vertLocToGloMap->GetLocalToGlobalMap();
// number of Dir coeffs in from full problem
int nDirFull = m_locToGloMap->GetNumGlobalDirBndCoeffs();
Array<OneD,NekDouble> In(nglo,0.0);
Array<OneD,NekDouble> Out(nglo,0.0);
// Gather rhs
for(i = 0; i < nloc; ++i)
{
val = LocToGloBnd[i];
if(val >= nDirFull)
{
In[LocToGlo[i]] = pInput[val-nDirFull];
}
}
// Do solve without enforcing any boundary conditions.
m_vertLinsys->SolveLinearSystem(m_vertLocToGloMap->GetNumLocalCoeffs(),
In,Out,m_vertLocToGloMap);
if(pNonVertOutput != NullNekDouble1DArray)
{
ASSERTL1(pNonVertOutput.num_elements() >= pOutput.num_elements(),"Non Vert output is not of sufficient length");
Vmath::Vcopy(pOutput.num_elements(),pNonVertOutput,1,pOutput,1);
}
else
{
//Copy input to output as a unit preconditioner on
//any other value
Vmath::Vcopy(pInput.num_elements(),pInput,1,pOutput,1);
}
// Scatter back soln from linear solve
for(i = 0; i < nloc; ++i)
{
val = LocToGloBnd[i];
if(val >= nDirFull)
{
pOutput[val-nDirFull] = Out[LocToGlo[i]];
}
}
}
break;
default:
ASSERTL0(0,"Unsupported solver type");
break;
}
}
}
}
///////////////////////////////////////////////////////////////////////////////
//
// File PreconditionerPETScLinear.h
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//