Commit 15c51a34 authored by Dave Moxey's avatar Dave Moxey

Merge branch 'master' into fix/CFS-mu-pu

parents fbfc62b3 36763032
......@@ -5,6 +5,9 @@
SET(MemoryManagerSources
MemoryManager.cpp)
SET(PartitionAnalyseSources
PartitionAnalyse.cpp)
SET(VecMatSources
VecMatExample.cpp)
......@@ -37,6 +40,9 @@ SET(TimeIntegrationDemoSources
ADD_NEKTAR_EXECUTABLE(Matrix demos VecMatSources)
SET_LAPACK_LINK_LIBRARIES(Matrix)
ADD_NEKTAR_EXECUTABLE(PartitionAnalyse demos PartitionAnalyseSources)
SET_LAPACK_LINK_LIBRARIES(PartitionAnalyse)
ADD_NEKTAR_EXECUTABLE(FoundationDemo demos FoundationSources )
SET_LAPACK_LINK_LIBRARIES(FoundationDemo)
......
///////////////////////////////////////////////////////////////////////////////
//
// File: PartitionAnalyse.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: Small utility to export histogram of partition sizes.
//
///////////////////////////////////////////////////////////////////////////////
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <LibUtilities/BasicUtils/MeshPartition.h>
#include <LibUtilities/Communication/CommSerial.h>
#include <iostream>
using namespace std;
using namespace Nektar;
using namespace Nektar::LibUtilities;
class FauxComm : public CommSerial
{
public:
FauxComm(int argc, char* argv[], int size)
: CommSerial(argc, argv)
{
m_size = size;
m_type = "Faux parallel";
}
FauxComm(int size) : CommSerial(0, NULL)
{
m_size = size;
m_type = "Faux parallel";
}
virtual ~FauxComm() {}
void v_SplitComm(int pRows, int pColumns)
{
m_commRow = boost::shared_ptr<FauxComm>(new FauxComm(pColumns));
m_commColumn = boost::shared_ptr<FauxComm>(new FauxComm(pRows));
}
};
int main(int argc, char *argv[])
{
if (argc < 3)
{
cerr << "Usage: PartitionAnalyse <nproc> <xml file1> [xml file 2..n]"
<< endl;
return 1;
}
int nParts = atoi(argv[1]);
vector<string> filenames(argv + 2, argv + argc);
CommSharedPtr vComm = boost::shared_ptr<FauxComm>(
new FauxComm(argc, argv, nParts));
char **new_argv = new char*[argc];
new_argv[0] = strdup("PartitionAnalyse");
new_argv[1] = strdup("--part-info");
for (int i = 0; i < argc-2; ++i)
{
new_argv[i+2] = strdup(filenames[i].c_str());
}
LibUtilities::SessionReaderSharedPtr vSession
= LibUtilities::SessionReader::CreateInstance(
argc, new_argv, filenames, vComm);
return 0;
}
......@@ -468,6 +468,105 @@ namespace Nektar
ParseUtils::GenerateSeqVector(vSeqStr.c_str(), m_domain);
}
void MeshPartition::PrintPartInfo(std::ostream &out)
{
int nElmt = boost::num_vertices(m_mesh);
int nPart = m_comm->GetRowComm()->GetSize();
out << "# Partition information:" << std::endl;
out << "# No. elements : " << nElmt << std::endl;
out << "# No. partitions: " << nPart << std::endl;
BoostVertexIterator vertit, vertit_end;
std::vector<int> partElmtCount(nPart, 0);
std::vector<int> partLocCount (nPart, 0);
std::vector<int> partBndCount (nPart, 0);
std::map<int, int> elmtSizes;
std::map<int, int> elmtBndSizes;
for (unsigned int i = 0; i < m_domain.size(); ++i)
{
int cId = m_domain[i];
NummodesPerField npf = m_expansions[cId];
for (NummodesPerField::iterator it = npf.begin(); it != npf.end(); ++it)
{
ASSERTL0(it->second.size() == m_dim,
" Number of directional" \
" modes in expansion spec for composite id = " +
boost::lexical_cast<std::string>(cId) +
" and field " +
boost::lexical_cast<std::string>(it->first) +
" does not correspond to mesh dimension");
int na = it->second[0];
int nb = it->second[1];
int nc = 0;
if (m_dim == 3)
{
nc = it->second[2];
}
int weight = 0;
int bndWeight = 0;
switch (m_meshComposites[cId].type)
{
case 'A':
weight = StdTetData::getNumberOfCoefficients(na, nb, nc);
bndWeight = StdTetData::getNumberOfBndCoefficients(na, nb, nc);
break;
case 'R':
weight = StdPrismData::getNumberOfCoefficients(na, nb, nc);
bndWeight = StdPrismData::getNumberOfBndCoefficients(na, nb, nc);
break;
case 'H':
weight = StdHexData::getNumberOfCoefficients(na, nb, nc);
bndWeight = StdHexData::getNumberOfBndCoefficients(na, nb, nc);
break;
case 'P':
weight = StdPyrData::getNumberOfCoefficients(na, nb, nc);
bndWeight = StdPyrData::getNumberOfBndCoefficients(na, nb, nc);
break;
case 'Q':
weight = StdQuadData::getNumberOfCoefficients(na, nb);
bndWeight = StdQuadData::getNumberOfBndCoefficients(na, nb);
break;
case 'T':
weight = StdTriData::getNumberOfCoefficients(na, nb);
bndWeight = StdTriData::getNumberOfBndCoefficients(na, nb);
break;
case 'S':
weight = StdSegData::getNumberOfCoefficients(na);
bndWeight = StdSegData::getNumberOfBndCoefficients(na);
break;
default:
break;
}
for (unsigned int j = 0; j < m_meshComposites[cId].list.size(); ++j)
{
int elid = m_meshComposites[cId].list[j];
elmtSizes[elid] = weight;
elmtBndSizes[elid] = bndWeight;
}
}
}
for (boost::tie(vertit, vertit_end) = boost::vertices(m_mesh);
vertit != vertit_end; ++vertit)
{
int partId = m_mesh[*vertit].partition;
partElmtCount[partId]++;
partLocCount [partId] += elmtSizes[m_mesh[*vertit].id];
partBndCount [partId] += elmtBndSizes[m_mesh[*vertit].id];
}
for (int i = 0; i < nPart; ++i)
{
out << i << " " << partElmtCount[i] << " " << partLocCount[i] << " " << partBndCount[i] << std::endl;
}
}
void MeshPartition::ReadConditions(const SessionReaderSharedPtr& pSession)
{
......@@ -540,7 +639,7 @@ namespace Nektar
*/
void MeshPartition::WeightElements()
{
std::vector<unsigned int> weight(2*m_numFields, 1);
std::vector<unsigned int> weight(m_numFields, 1);
for (int i = 0; i < m_meshElements.size(); ++i)
{
m_vertWeights.push_back( weight );
......@@ -562,33 +661,48 @@ namespace Nektar
" does not correspond to mesh dimension");
int na = it->second[0];
int nb = it->second[1];
int nb = 0;
int nc = 0;
if (m_dim >= 2)
{
nb = it->second[1];
}
if (m_dim == 3)
{
nc = it->second[2];
}
int weight = 0;
int bndWeight = 0;
switch (m_meshComposites[cId].type)
{
case 'A':
weight = StdTetData::getNumberOfCoefficients(na, nb, nc);
weight = StdTetData::getNumberOfCoefficients(na, nb, nc);
bndWeight = StdTetData::getNumberOfBndCoefficients(na, nb, nc);
break;
case 'R':
weight = StdPrismData::getNumberOfCoefficients(na, nb, nc);
case 'R':
weight = StdPrismData::getNumberOfCoefficients(na, nb, nc);
bndWeight = StdPrismData::getNumberOfBndCoefficients(na, nb, nc);
break;
case 'H':
weight = StdHexData::getNumberOfCoefficients(na, nb, nc);
case 'H':
weight = StdHexData::getNumberOfCoefficients(na, nb, nc);
bndWeight = StdHexData::getNumberOfBndCoefficients(na, nb, nc);
break;
case 'P':
weight = StdPyrData::getNumberOfCoefficients(na, nb, nc);
case 'P':
weight = StdPyrData::getNumberOfCoefficients(na, nb, nc);
bndWeight = StdPyrData::getNumberOfBndCoefficients(na, nb, nc);
break;
case 'Q':
weight = StdQuadData::getNumberOfCoefficients(na, nb);
case 'Q':
weight = StdQuadData::getNumberOfCoefficients(na, nb);
bndWeight = StdQuadData::getNumberOfBndCoefficients(na, nb);
break;
case 'T':
weight = StdTriData::getNumberOfCoefficients(na, nb);
case 'T':
weight = StdTriData::getNumberOfCoefficients(na, nb);
bndWeight = StdTriData::getNumberOfBndCoefficients(na, nb);
break;
case 'S':
weight = StdSegData::getNumberOfCoefficients(na);
bndWeight = StdSegData::getNumberOfBndCoefficients(na);
break;
default:
break;
......@@ -597,8 +711,8 @@ namespace Nektar
for (unsigned int j = 0; j < m_meshComposites[cId].list.size(); ++j)
{
int elmtId = m_meshComposites[cId].list[j];
m_vertWeights[elmtId][ m_fieldNameToId[ it->first ] * 2 + 0 ] = weight;
m_vertWeights[elmtId][ m_fieldNameToId[ it->first ] * 2 + 1 ] = weight*weight;
m_vertWeights[elmtId][ m_fieldNameToId[it->first]] = bndWeight;
//m_vertWeights[elmtId][ m_fieldNameToId[ it->first ] * 2 + 1 ] = weight*weight;
}
}
} // for i
......@@ -656,8 +770,7 @@ namespace Nektar
{
int acnt = 0;
int vcnt = 0;
int nWeight = m_weightingRequired ? 2*nGraphVerts*m_numFields
: nGraphVerts;
int nWeight = nGraphVerts;
BoostAdjacencyIterator adjvertit, adjvertit_end;
Array<OneD, int> xadj(nGraphVerts+1,0);
Array<OneD, int> adjncy(2*nGraphEdges);
......@@ -678,11 +791,7 @@ namespace Nektar
if (m_weightingRequired)
{
// populate vertex multi-weights
for (i = 0; i < 2*m_numFields; i++)
{
vwgt[pGraph[*vertit].id * m_numFields + i] = pGraph[*vertit].weight[i];
}
vwgt[pGraph[*vertit].id ] = pGraph[*vertit].weight[0];
}
else
{
......@@ -702,7 +811,7 @@ namespace Nektar
if(m_comm->GetColumnComm()->GetRank() == 0)
{
// Attempt partitioning using METIS.
int ncon = m_weightingRequired ? 2*m_numFields : 1;
int ncon = 1;
Metis::PartGraphVKway(nGraphVerts, ncon, xadj, adjncy, vwgt, vsize, npart, vol, part);
// Check METIS produced a valid partition and fix if not.
CheckPartitions(part);
......
......@@ -64,6 +64,7 @@ namespace Nektar
LIB_UTILITIES_EXPORT void WriteAllPartitions(
SessionReaderSharedPtr& pSession);
LIB_UTILITIES_EXPORT void PrintPartInfo(std::ostream &out);
LIB_UTILITIES_EXPORT void GetCompositeOrdering(
CompositeOrdering &composites);
LIB_UTILITIES_EXPORT void GetBndRegionOrdering(
......
......@@ -222,6 +222,12 @@ namespace Nektar
else
{
m_comm = pComm;
if (m_comm->GetSize() > 1)
{
m_solverInfoDefaults["GLOBALSYSSOLN"] =
"IterativeStaticCond";
}
}
// If running in parallel change the default global sys solution
......@@ -289,7 +295,7 @@ namespace Nektar
"number of procs in Y-dir")
("npz", po::value<int>(),
"number of procs in Z-dir")
("part-info", "Output partition information")
;
CmdLineArgMap::const_iterator cmdIt;
......@@ -1367,6 +1373,11 @@ namespace Nektar
vPartitioner->WriteAllPartitions(vSession);
vPartitioner->GetCompositeOrdering(m_compOrder);
vPartitioner->GetBndRegionOrdering(m_bndRegOrder);
if (DefinesCmdLineArgument("part-info"))
{
vPartitioner->PrintPartInfo(std::cout);
}
}
}
else
......@@ -1383,6 +1394,11 @@ namespace Nektar
vPartitioner->WriteLocalPartition(vSession);
vPartitioner->GetCompositeOrdering(m_compOrder);
vPartitioner->GetBndRegionOrdering(m_bndRegOrder);
if (DefinesCmdLineArgument("part-info"))
{
vPartitioner->PrintPartInfo(std::cout);
}
}
m_comm->Block();
......
......@@ -92,23 +92,56 @@ namespace Nektar
};
namespace StdSegData
{
inline int getNumberOfCoefficients(int Na)
{
return Na;
}
inline int getNumberOfBndCoefficients(int Na)
{
return 2;
}
}
// Dimensions of coefficients for each space
namespace StdTriData
{
inline int getNumberOfCoefficients(int Na, int Nb)
{
ASSERTL0(Na <= Nb, "order in 'a' direction is higher "
ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL1(Na <= Nb, "order in 'a' direction is higher "
"than order in 'b' direction");
return Na*(Na+1)/2 + Na*(Nb-Na);
}
inline int getNumberOfBndCoefficients(int Na, int Nb)
{
ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL1(Na <= Nb, "order in 'a' direction is higher "
"than order in 'b' direction");
return (Na-1) + 2*(Nb-1);
}
}
namespace StdQuadData
{
inline int getNumberOfCoefficients(int Na, int Nb)
{
ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
return Na*Nb;
}
inline int getNumberOfBndCoefficients(int Na, int Nb)
{
ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
return 2*(Na-1) + 2*(Nb-1);
}
}
......@@ -117,8 +150,20 @@ namespace Nektar
{
inline int getNumberOfCoefficients( int Na, int Nb, int Nc )
{
ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL2(Nc > 1, "Order in 'c' direction must be > 1.");
return Na*Nb*Nc;
}
inline int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
{
ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL2(Nc > 1, "Order in 'c' direction must be > 1.");
return 2*Na*Nb + 2*Na*Nc + 2*Nb*Nc
- 4*(Na + Nb + Nc) + 8;
}
}
namespace StdTetData
......@@ -140,6 +185,13 @@ namespace Nektar
*/
inline int getNumberOfCoefficients(int Na, int Nb, int Nc)
{
ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL2(Nc > 1, "Order in 'c' direction must be > 1.");
ASSERTL1(Na <= Nb, "order in 'a' direction is higher "
"than order in 'b' direction");
ASSERTL1(Nb <= Nc, "order in 'b' direction is higher "
"than order in 'c' direction");
int nCoef = 0;
for (int a = 0; a < Na; ++a)
{
......@@ -153,6 +205,25 @@ namespace Nektar
}
return nCoef;
}
inline int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
{
ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL2(Nc > 1, "Order in 'c' direction must be > 1.");
ASSERTL1(Na <= Nb, "order in 'a' direction is higher "
"than order in 'b' direction");
ASSERTL1(Nb <= Nc, "order in 'b' direction is higher "
"than order in 'c' direction");
int nCoef = Na*(Na+1)/2 + (Nb-Na)*Na // base
+ Na*(Na+1)/2 + (Nc-Na)*Na // front
+ 2*(Nb*(Nb+1)/2 + (Nc-Nb)*Nb)// 2 other sides
- Na - 2*Nb - 3*Nc // less edges
+ 4; // plus vertices
return nCoef;
}
}
......@@ -160,6 +231,13 @@ namespace Nektar
{
inline int getNumberOfCoefficients(int Na, int Nb, int Nc)
{
ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
"than order in 'c' direction.");
ASSERTL1(Nb <= Nc, "Order in 'b' direction is higher "
"than order in 'c' direction.");
int nCoef = 0;
for (int c = 0; c < Nc; ++c)
{
......@@ -171,31 +249,54 @@ namespace Nektar
}
}
}
/*
for (int a = 0; a < Na; ++a)
{
for (int b = 0; b < Nb; ++b)
{
for (int c = 0; c < Nc - a - b; ++c)
{
++nCoef;
}
}
}
*/
//std::cout << "Na = " << Na << " Nb = " << Nb << " Nc = " << Nc << " nCoef = " << nCoef << std::endl;
return nCoef;
}
inline int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
{
ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
"than order in 'c' direction.");
ASSERTL1(Nb <= Nc, "Order in 'b' direction is higher "
"than order in 'c' direction.");
return Na*Nb // base
+ 2*(Na*(Na+1)/2 + (Nc-Na)*Na) // front and back
+ 2*(Nb*(Nb+1)/2 + (Nc-Nb)*Nb) // sides
- 2*Na - 2*Nb - 4*Nc // less edges
+ 5; // plus vertices
}
}
namespace StdPrismData
{
inline int getNumberOfCoefficients( int Na, int Nb, int Nc )
{
ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
"than order in 'c' direction.");
return Nb*StdTriData::getNumberOfCoefficients(Na,Nc);
}
}
inline int getNumberOfBndCoefficients( int Na, int Nb, int Nc)
{
ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
"than order in 'c' direction.");
return Na*Nb + 2*Nb*Nc // rect faces
+ 2*( Na*(Na+1)/2 + (Nc-Na)*Na ) // tri faces
- 2*Na - 3*Nb - 4*Nc // less edges
+ 6; // plus vertices
}
}
inline int GetNumberOfCoefficients(ShapeType shape, std::vector<unsigned int> &modes, int offset)
{
......
......@@ -1220,60 +1220,92 @@ namespace Nektar
Array<OneD, NekDouble> &locCoords,
NekDouble tol)
{
static int start = 0;
NekDouble resid, min_resid = NekConstants::kNekMinResidInit;
int min_elmt;
Array<OneD, NekDouble> min_locCoords(locCoords.num_elements());
NekDouble resid;
// start search at previous element or 0
for (int i = start; i < (*m_exp).size(); ++i)
if (GetNumElmts() == 0)
{
if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords, locCoords,
tol, resid))
return -1;
}
// Manifold case (point may match multiple elements)
if (GetExp(0)->GetCoordim() > GetExp(0)->GetShapeDimension())
{
std::vector<std::pair<int,NekDouble> > elmtIdDist;
SpatialDomains::PointGeomSharedPtr v;
SpatialDomains::PointGeom w;