Commit 58010c0f authored by Pavel Burovskiy's avatar Pavel Burovskiy
Browse files

Merge branch 'master' into fix/blas-calls-in-SL-staticcond

parents 64418b55 aaecbf1c
......@@ -160,6 +160,7 @@ INCLUDE (ThirdPartyFFTW)
INCLUDE (ThirdPartyArpack)
INCLUDE (ThirdPartyMPI)
INCLUDE (ThirdPartyVTK)
INCLUDE (ThirdPartyQT4)
IF( NEKTAR_USE_MKL )
INCLUDE (FindMKL)
......
SET(METIS_SEARCH_PATHS
${CMAKE_SOURCE_DIR}/ThirdParty/modmetis-4.0/
${CMAKE_SOURCE_DIR}/ThirdParty/modmetis-4.0/build/
${CMAKE_SOURCE_DIR}/../ThirdParty/modmetis-4.0/
${CMAKE_SOURCE_DIR}/../ThirdParty/modmetis-4.0/build
${CMAKE_SOURCE_DIR}/ThirdParty/metis-5.0.2/
${CMAKE_SOURCE_DIR}/ThirdParty/metis-5.0.2/build/
${CMAKE_SOURCE_DIR}/../ThirdParty/metis-5.0.2/
${CMAKE_SOURCE_DIR}/../ThirdParty/metis-5.0.2/build
${CMAKE_SOURCE_DIR}/ThirdParty/dist/lib
${CMAKE_SOURCE_DIR}/../ThirdParty/dist/lib)
FIND_LIBRARY(METIS_LIB NAMES modmetis PATHS ${METIS_SEARCH_PATHS})
FIND_LIBRARY(METIS_LIB NAMES metis PATHS ${METIS_SEARCH_PATHS})
SET(METIS_FOUND FALSE)
IF (METIS_LIB)
SET(METIS_FOUND TRUE)
GET_FILENAME_COMPONENT(METIS_PATH ${METIS_LIB} PATH)
INCLUDE_DIRECTORIES(${METIS_PATH}/../include)
MARK_AS_ADVANCED(METIS_PATH)
MARK_AS_ADVANCED(METIS_LIB)
ENDIF (METIS_LIB)
......
......@@ -221,7 +221,7 @@ MACRO(ADD_NEKTAR_LIBRARY name component type)
# NIST Sparse BLAS only static, so link into Nektar libraries directly.
TARGET_LINK_LIBRARIES( ${name} ${NIST_SPARSE_BLAS} ${METIS_LIB})
ADD_DEPENDENCIES(${name} spblastk0.9b modmetis-4.0 boost tinyxml zlib)
ADD_DEPENDENCIES(${name} spblastk0.9b modmetis-5.0.2 boost tinyxml zlib)
SET_PROPERTY(TARGET ${name} PROPERTY FOLDER ${component})
IF (NEKTAR_USE_MPI)
TARGET_LINK_LIBRARIES( ${name} ${GSMPI_LIBRARY} ${XXT_LIBRARY})
......
......@@ -4,17 +4,18 @@ SET(THIRDPARTY_BUILD_METIS ON CACHE BOOL
IF (THIRDPARTY_BUILD_METIS)
INCLUDE(ExternalProject)
EXTERNALPROJECT_ADD(
modmetis-4.0
modmetis-5.0.2
PREFIX ${TPSRC}
URL ${TPURL}/modmetis-4.0.tar.bz2
URL_MD5 "4267355e04dcc2ef36f9889a98c923a5"
URL ${TPURL}/modmetis-5.0.2.tar.bz2
URL_MD5 "ffbdc6a50283934389a0b3b0c32b62c0"
DOWNLOAD_DIR ${TPSRC}
CONFIGURE_COMMAND ${CMAKE_COMMAND} -DCMAKE_INSTALL_PREFIX:PATH=${TPSRC}/dist ${TPSRC}/src/modmetis-4.0
CONFIGURE_COMMAND ${CMAKE_COMMAND} -DCMAKE_INSTALL_PREFIX:PATH=${TPSRC}/dist -DCMAKE_C_FLAGS:STRING=-fPIC -DGKLIB_PATH:PATH=${TPSRC}/src/modmetis-5.0.2/GKlib ${TPSRC}/src/modmetis-5.0.2
)
SET(METIS_LIB modmetis CACHE FILEPATH
SET(METIS_LIB metis CACHE FILEPATH
"METIS library" FORCE)
MARK_AS_ADVANCED(METIS_LIB)
LINK_DIRECTORIES(${TPSRC}/dist/lib)
INCLUDE_DIRECTORIES(${TPSRC}/dist/include)
ELSE (THIRDPARTY_BUILD_METIS)
INCLUDE (FindMetis)
ENDIF (THIRDPARTY_BUILD_METIS)
......
# QT4
OPTION(NEKTAR_USE_QT4
"Use QT4 for graphical tools and utilities." OFF)
IF( NEKTAR_USE_QT4 )
FIND_PACKAGE(Qt4 COMPONENTS QtCore QtGui REQUIRED)
INCLUDE(${QT_USE_FILE})
ADD_DEFINITIONS(${QT_DEFINITIONS})
ENDIF( NEKTAR_USE_QT4 )
......@@ -8,6 +8,7 @@ library/CPackRPM.cmake
library/Demos/StdRegions/ExtraDemos
library/Demos/MultiRegions/ExtraDemos
solvers/ADR2DManifoldSolver
solvers/CardiacEPSolver/Utilities
solvers/FitzHughNagumoSolver
solvers/ImageWarpingSolver
solvers/PulseWaveSolver
......
......@@ -341,7 +341,13 @@ namespace Nektar
try
{
// Attempt partitioning using METIS.
Metis::PartGraphVKway(nGraphVerts, xadj, adjncy, vwgt, vsize, npart, vol, part);
// Check METIS produced a valid partition and fix if not.
CheckPartitions(part);
// Distribute partitioning to all processes.
for (i = 1; i < m_comm->GetSize(); ++i)
{
m_comm->Send(i, part);
......@@ -376,6 +382,39 @@ namespace Nektar
}
}
void MeshPartition::CheckPartitions(Array<OneD, int> &pPart)
{
unsigned int i = 0;
unsigned int cnt = 0;
const unsigned int npart = m_comm->GetSize();
bool valid = true;
// Check that every process has at least one element assigned
for (i = 0; i < npart; ++i)
{
cnt = std::count(pPart.begin(), pPart.end(), i);
if (cnt == 0)
{
valid = false;
}
}
// If METIS produced an invalid partition, repartition naively.
// Elements are assigned to processes in a round-robin fashion.
// It is assumed that METIS failure only occurs when the number of
// elements is approx. the number of processes, so this approach
// should not be too inefficient communication-wise.
if (!valid)
{
for (i = 0; i < pPart.num_elements(); ++i)
{
pPart[i] = i % npart;
}
}
}
void MeshPartition::OutputPartition(
LibUtilities::SessionReaderSharedPtr& pSession,
BoostSubGraph& pGraph,
......
......@@ -180,6 +180,7 @@ namespace Nektar
void PartitionGraph(BoostSubGraph& pGraph,
BoostSubGraph& pLocalPartition);
void OutputPartition(SessionReaderSharedPtr& pSession, BoostSubGraph& pGraph, TiXmlElement* pGeometry);
void CheckPartitions(Array<OneD, int> &pPart);
};
typedef boost::shared_ptr<MeshPartition> MeshPartitionSharedPtr;
......
......@@ -39,99 +39,63 @@
#include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include "metis.h"
namespace Metis
{
extern "C"
{
// -- Sparse MAtrix Reordering (equivalent to onmetis)
void METIS_NodeND (int *nVerts, int *xadj, int *adjncy, int *numflag, int *options,
int *perm, int *iperm);
void AS_METIS_NodeND(int *nVerts, int *xadj, int *adjncy, int *numflag, int *options,
int *perm, int *iperm, int *map, int *mdswitch);
void METIS_PartMeshNodal(int *nElmts, int *nVerts, int *elType, int *numflag,
int *nparts, int *edgecut, int *epart, int *npart);
void METIS_PartGraphVKway(int *nVerts, int *xadj, int *adjcy, int *vertWgt, int *vertSize,
int *wgtFlag, int *numflag, int *nparts, int *options, int *volume, int *part);
}
//#ifdef NEKTAR_USING_METIS
inline static void onmetis(int *nVerts, int *xadj, int *adjncy, int *numflag, int *options,
int *perm, int *iperm)
{
METIS_NodeND(nVerts,xadj,adjncy,numflag,options,perm,iperm);
}
inline static void as_onmetis(int *nVerts, int *xadj, int *adjncy, int *numflag, int *options,
int *perm, int *iperm, int *map, int *mdswitch)
{
AS_METIS_NodeND(nVerts,xadj,adjncy,numflag,options,perm,iperm,map,mdswitch) ;
void AS_METIS_NodeND(int *nVerts, int *xadj, int *adjncy, int *vwgt,
int *options, int *perm, int *iperm, int *map,
int *mdswitch);
}
inline static void onmetis(int nVerts, Nektar::Array<Nektar::OneD, int> xadj, Nektar::Array<Nektar::OneD, int> adjncy,
Nektar::Array<Nektar::OneD, int> perm, Nektar::Array<Nektar::OneD, int> iperm)
inline static void as_onmetis(
int nVerts,
Nektar::Array<Nektar::OneD, int> xadj,
Nektar::Array<Nektar::OneD, int> adjncy,
Nektar::Array<Nektar::OneD, int> perm,
Nektar::Array<Nektar::OneD, int> iperm,
Nektar::Array<Nektar::OneD, int> map,
int mdswitch = 1)
{
ASSERTL1(xadj.num_elements() == nVerts+1,"Array xadj out of bounds");
ASSERTL1(perm.num_elements() == nVerts,"Array perm out of bounds");
ASSERTL1(iperm.num_elements() == nVerts,"Array iperm out of bounds");
int numflag = 0;
int options[8];
options[0]=0;
METIS_NodeND(&nVerts,&xadj[0],&adjncy[0],&numflag,options,&perm[0],&iperm[0]);
AS_METIS_NodeND(&nVerts, &xadj[0], &adjncy[0], NULL, NULL, &perm[0],
&iperm[0], &map[0], &mdswitch);
}
inline static void as_onmetis(int nVerts, Nektar::Array<Nektar::OneD, int> xadj, Nektar::Array<Nektar::OneD, int> adjncy,
Nektar::Array<Nektar::OneD, int> perm, Nektar::Array<Nektar::OneD, int> iperm, Nektar::Array<Nektar::OneD, int> map,
int mdswitch)
{
ASSERTL1(xadj.num_elements() == nVerts+1,"Array xadj out of bounds");
ASSERTL1(perm.num_elements() == nVerts,"Array perm out of bounds");
ASSERTL1(iperm.num_elements() == nVerts,"Array iperm out of bounds");
int numflag = 0;
int options[8];
options[0]=0;
AS_METIS_NodeND(&nVerts,&xadj[0],&adjncy[0],&numflag,options,&perm[0],&iperm[0],&map[0],&mdswitch);
}
inline static void MeshPartition(int nElmts, int nVerts, Nektar::Array<Nektar::OneD, int>& mesh, int type, int nparts,
Nektar::Array<Nektar::OneD, int>& edgePart, Nektar::Array<Nektar::OneD, int>& nodePart)
{
int numflag = 0;
METIS_PartMeshNodal(&nElmts, &nVerts, &mesh[0], &type, &numflag, &nparts, &edgePart[0], &nodePart[0]);
}
inline static void PartGraphVKway( int& nVerts,
Nektar::Array<Nektar::OneD, int>& xadj,
Nektar::Array<Nektar::OneD, int>& adjcy,
Nektar::Array<Nektar::OneD, int>& vertWgt,
Nektar::Array<Nektar::OneD, int>& vertSize,
int& nparts,
int& volume,
Nektar::Array<Nektar::OneD, int>& part)
inline static void PartGraphVKway(
int& nVerts,
Nektar::Array<Nektar::OneD, int>& xadj,
Nektar::Array<Nektar::OneD, int>& adjcy,
Nektar::Array<Nektar::OneD, int>& vertWgt,
Nektar::Array<Nektar::OneD, int>& vertSize,
int& nparts,
int& volume,
Nektar::Array<Nektar::OneD, int>& part)
{
int wgtflag = 0;
int *wgts = 0;
int *sizes = 0;
int *vwgt = 0;
int *vsize = 0;
if (vertWgt.num_elements() > 0)
{
wgtflag += 1;
wgts = &vertWgt[0];
vwgt = &vertWgt[0];
}
if (vertSize.num_elements() > 0)
{
wgtflag += 2;
sizes = &vertSize[0];
vsize = &vertSize[0];
}
int numflag = 0;
int options[5];
options[0]=0;
METIS_PartGraphVKway(&nVerts, &xadj[0], &adjcy[0], wgts, sizes, &wgtflag,
&numflag, &nparts, options, &volume, &part[0]);
int ncon = 1;
int options[METIS_NOPTIONS];
METIS_SetDefaultOptions(options);
METIS_PartGraphKway(&nVerts, &ncon, &xadj[0], &adjcy[0], vwgt, vsize,
0, &nparts, 0, 0, options, &volume, &part[0]);
}
//#endif //NEKTAR_USING_METIS
}
#endif //NEKTAR_LIB_UTILITIES_BASICUTILS_METIS_HPP
......@@ -62,6 +62,32 @@ namespace Nektar
const Array<OneD, const Array<OneD, NekDouble> > normals
= GetEdgeNormal(edge);
if (m_requireNeg.size() == 0)
{
m_requireNeg.resize(GetNedges());
for (int i = 0; i < GetNedges(); ++i)
{
m_requireNeg[i] = false;
if (m_negatedNormals[i])
{
m_requireNeg[i] = true;
}
Expansion1DSharedPtr edgeExp = boost::dynamic_pointer_cast<
Expansion1D>(m_edgeExp[i].lock());
if (edgeExp->GetRightAdjacentElementExp())
{
if (edgeExp->GetRightAdjacentElementExp()->GetGeom2D()
->GetGlobalID() == GetGeom2D()->GetGlobalID())
{
m_requireNeg[i] = true;
}
}
}
}
// We allow the case of mixed polynomial order by supporting only
// those modes on the edge common to both adjoining elements. This
// is enforced here by taking the minimum size and padding with
......@@ -79,6 +105,7 @@ namespace Nektar
boost::dynamic_pointer_cast<
LocalRegions::Expansion1D>(EdgeExp);
if (m_negatedNormals[edge])
{
Vmath::Neg(nquad_e,EdgeExp->UpdatePhys(),1);
......@@ -104,6 +131,32 @@ namespace Nektar
{
int i;
if (m_requireNeg.size() == 0)
{
m_requireNeg.resize(GetNedges());
for (i = 0; i < GetNedges(); ++i)
{
m_requireNeg[i] = false;
if (m_negatedNormals[i])
{
m_requireNeg[i] = true;
}
Expansion1DSharedPtr edgeExp = boost::dynamic_pointer_cast<
Expansion1D>(m_edgeExp[i].lock());
if (edgeExp->GetRightAdjacentElementExp())
{
if (edgeExp->GetRightAdjacentElementExp()->GetGeom2D()
->GetGlobalID() == GetGeom2D()->GetGlobalID())
{
m_requireNeg[i] = true;
}
}
}
}
StdRegions::Orientation edgedir = GetEorient(edge);
unsigned short num_mod0 = EdgeExp->GetBasis(0)->GetNumModes();
......@@ -126,18 +179,10 @@ namespace Nektar
boost::dynamic_pointer_cast<
LocalRegions::Expansion1D>(EdgeExp);
if (m_negatedNormals[edge])
if (m_requireNeg[edge])
{
Vmath::Neg(n_coeffs,EdgeExp->UpdateCoeffs(),1);
}
else if (locExp->GetRightAdjacentElementEdge() != -1)
{
if (locExp->GetRightAdjacentElementExp()->GetGeom2D()->GetGlobalID()
== GetGeom2D()->GetGlobalID())
{
Vmath::Neg(n_coeffs,EdgeExp->UpdateCoeffs(),1);
}
}
Array<OneD, NekDouble> coeff(n_coeffs,0.0);
LibUtilities::BasisType btype = ((LibUtilities::BasisType) 1); //1-->Ortho_A
......@@ -162,18 +207,10 @@ namespace Nektar
boost::dynamic_pointer_cast<
LocalRegions::Expansion1D>(EdgeExp);
if (m_negatedNormals[edge])
if (m_requireNeg[edge])
{
Vmath::Neg(n_coeffs,EdgeExp->UpdateCoeffs(),1);
}
else if (locExp->GetRightAdjacentElementEdge() != -1)
{
if (locExp->GetRightAdjacentElementExp()->GetGeom2D()->GetGlobalID()
== GetGeom2D()->GetGlobalID())
{
Vmath::Neg(order_e,EdgeExp->UpdateCoeffs(),1);
}
}
}
// add data to outarray if forward edge normal is outwards
......
......@@ -134,6 +134,7 @@ namespace Nektar
private:
std::vector<ExpansionWeakPtr> m_edgeExp;
std::vector<bool> m_requireNeg;
Expansion3DWeakPtr m_elementLeft;
Expansion3DWeakPtr m_elementRight;
......
......@@ -42,7 +42,7 @@ namespace Nektar
{
namespace LocalRegions
{
Expansion3D::Expansion3D(){}
Expansion3D::Expansion3D() : m_requireNeg() {}
void Expansion3D::AddHDGHelmholtzTraceTerms(
const NekDouble tau,
......@@ -914,6 +914,40 @@ namespace Nektar
StdRegions::IndexMapValuesSharedPtr map =
StdExpansion::GetIndexMap(ikey);
/*
* Coming into this routine, the velocity V will have been
* multiplied by the trace normals to give the input vector Vn. By
* convention, these normals are inwards facing for elements which
* have FaceExp as their right-adjacent face. This conditional
* statement therefore determines whether the normals must be
* negated, since the integral being performed here requires an
* outwards facing normal.
*/
if (m_requireNeg.size() == 0)
{
m_requireNeg.resize(GetNfaces());
for (i = 0; i < GetNfaces(); ++i)
{
m_requireNeg[i] = false;
if (m_negatedNormals[i])
{
m_requireNeg[i] = true;
}
Expansion2DSharedPtr faceExp = m_faceExp[i].lock();
if (faceExp->GetRightAdjacentElementExp())
{
if (faceExp->GetRightAdjacentElementExp()->GetGeom3D()->GetGlobalID()
== GetGeom3D()->GetGlobalID())
{
m_requireNeg[i] = true;
}
}
}
}
int order_e = (*map).num_elements(); // Order of the element
int n_coeffs = FaceExp->GetCoeffs().num_elements(); // Order of the trace
......@@ -924,37 +958,21 @@ namespace Nektar
else
{
FaceExp->IProductWRTBase(Fn,FaceExp->UpdateCoeffs());
LocalRegions::Expansion2DSharedPtr locExp =
boost::dynamic_pointer_cast<
LocalRegions::Expansion2D>(FaceExp);
/*
* Coming into this routine, the velocity V will have been
* multiplied by the trace normals to give the input vector
* Vn. By convention, these normals are inwards facing for
* elements which have FaceExp as their right-adjacent face.
* This conditional statement therefore determines whether the
* normals must be negated, since the integral being performed
* here requires an outwards facing normal.
*/
if (m_negatedNormals[face])
{
Vmath::Neg(order_e,FaceExp->UpdateCoeffs(),1);
}
else if (locExp->GetRightAdjacentElementFace() != -1)
}
if (m_requireNeg[face])
{
for(i = 0; i < order_e; ++i)
{
if (locExp->GetRightAdjacentElementExp()->GetGeom3D()->GetGlobalID()
== GetGeom3D()->GetGlobalID())
{
Vmath::Neg(order_e,FaceExp->UpdateCoeffs(),1);
}
outarray[(*map)[i].index] -= (*map)[i].sign*FaceExp->GetCoeff(i);
}
}
for(i = 0; i < order_e; ++i)
else
{
outarray[(*map)[i].index] += (*map)[i].sign*FaceExp->GetCoeff(i);
for(i = 0; i < order_e; ++i)
{
outarray[(*map)[i].index] += (*map)[i].sign*FaceExp->GetCoeff(i);
}
}
}
......
......@@ -106,6 +106,7 @@ namespace Nektar
// Only use this class for member functions
std::vector<Expansion2DWeakPtr> m_faceExp;
std::vector<bool> m_requireNeg;
};
// type defines for use of PrismExp in a boost vector
......
......@@ -368,24 +368,22 @@ namespace Nektar
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
{
const int nqtot = GetTotPoints();
Array<OneD, const NekDouble> jac = m_metricinfo->GetJac();
Array<OneD, NekDouble> tmp(nqtot);
// multiply inarray with Jacobian
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
Vmath::Vmul(nqtot, &jac[0], 1, (NekDouble*)&inarray[0], 1,
&tmp[0], 1);
}
else
{
Vmath::Smul(nqtot, jac[0], (NekDouble*)&inarray[0], 1,
&tmp[0], 1);
}
StdHexExp::v_IProductWRTBase(tmp, outarray);
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
int order0 = m_base[0]->GetNumModes();
int order1 = m_base[1]->GetNumModes();
Array<OneD, NekDouble> tmp(inarray.num_elements());
Array<OneD, NekDouble> wsp(nquad0*nquad1*(nquad2+order0) +
order0*order1*nquad2);
MultiplyByQuadratureMetric(inarray, tmp);
IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
m_base[1]->GetBdata(),
m_base[2]->GetBdata(),
tmp,outarray,wsp,
true,true,true);
}
void HexExp::v_IProductWRTDerivBase(
......@@ -433,50 +431,49 @@ namespace Nektar
const Array<TwoD, const NekDouble>& gmat = m_metricinfo->GetGmat();
Array<OneD, NekDouble> alloc(3*nqtot + 2*m_ncoeffs +
Array<OneD, NekDouble> alloc(4*nqtot + 2*m_ncoeffs +
nmodes0*nquad2*(nquad1+nmodes1));
Array<OneD, NekDouble> tmp1(alloc); // Dir1 metric
Array<OneD, NekDouble> tmp2(alloc + nqtot); // Dir2 metric
Array<OneD, NekDouble> tmp3(alloc + 2*nqtot); // Dir3 metric
Array<OneD, NekDouble> tmp4(alloc + 3*nqtot); // Dir1 iprod
Array<OneD, NekDouble> tmp5(tmp4 + m_ncoeffs); // Dir2 iprod
Array<OneD, NekDouble> wsp (tmp4 + 2*m_ncoeffs); // Wsp
Array<OneD, NekDouble> tmp1(alloc); // Quad metric
Array<OneD, NekDouble> tmp2(alloc + nqtot); // Dir1 metric
Array<OneD, NekDouble> tmp3(alloc + 2*nqtot); // Dir2 metric
Array<OneD, NekDouble> tmp4(alloc + 3*nqtot); // Dir3 metric
Array<OneD, NekDouble> tmp5(alloc + 4*nqtot); // Dir1 iprod
Array<OneD, NekDouble> tmp6(tmp5 + m_ncoeffs); // Dir2 iprod
Array<OneD, NekDouble> wsp (tmp5 + 2*m_ncoeffs); // Wsp
MultiplyByQuadratureMetric(inarray, tmp1);
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
Vmath::Vmul(nqtot,&gmat[3*dir][0], 1,inarray.get(),1,tmp1.get(),1);
Vmath::Vmul(nqtot,&gmat[3*dir+1][0],1,inarray.get(),1,tmp2.get(),1);
Vmath::Vmul(nqtot,&gmat[3*dir+2][0],1,inarray.get(),1,tmp3.get(),1);