parallel programming

This commit is contained in:
Joachim Schoeberl 2011-02-14 12:27:18 +00:00
parent 789b56179e
commit 846542c45c
17 changed files with 253 additions and 128 deletions

View File

@ -2,8 +2,8 @@ ACLOCAL_AMFLAGS = -I m4
METASOURCES = AUTO
SUBDIRS = libsrc ng nglib tutorials doc windows
SUBDIRS = libsrc ng tutorials doc windows
# nglib
# TESTS = ng/netgen -batchmode

View File

@ -7,7 +7,7 @@ AC_PREFIX_DEFAULT(["/opt/netgen"])
# Tcl/Tk configuration:
TEA_INIT([3.6])
TEA_INIT([3.9])
TEA_PATH_TCLCONFIG
TEA_LOAD_TCLCONFIG
TEA_PATH_TKCONFIG
@ -89,11 +89,15 @@ AC_ARG_ENABLE([nglib],
AC_ARG_ENABLE([parallel],
[AC_HELP_STRING([--enable-parallel],[enable mpi parallelization])],
[AC_SUBST([MPI_INCLUDES], "-I/opt/mpich/include -DPARALLEL -I/usr/share/metis/Lib -DMETIS")
AC_SUBST([MPI_LIBS], "-lmetis -L/opt/mpich/ch-p4/lib -lmpich")
[AC_SUBST([MPI_INCLUDES], "-I/usr/lib64/mpi/gcc/openmpi/include -DPARALLEL -I/usr/include/metis -DMETIS")
AC_SUBST([MPI_LIBS], "-L/usr/lib64/mpi/gcc/openmpi/lib64 -lmpi_cxx -lmetis")
]
)
# [AC_SUBST([MPI_INCLUDES], "-I/opt/mpich/include -DPARALLEL -I/usr/include/metis -DMETIS")
# AC_SUBST([MPI_LIBS], " -L/opt/mpich/ch-p4/lib64 -lmpich -lmetis")
# [AC_SUBST([MPI_INCLUDES], "-DPARALLEL -DNO_PARALLEL_THREADS -I/opt/mpich/include")
# AC_SUBST([MPI_LIBS], "")
# -lmetis
# [AC_SUBST([MPI_INCLUDES], "-I/opt/mpich/include -DPARALLEL -I/home/joachim/download/metis-4.0/Lib -DMETIS")
# AC_SUBST([MPI_LIBS], "-L/home/joachim/download/metis-4.0 -lmetis -L/opt/mpich/ch-p4/lib64 -lmpich")

View File

@ -41,7 +41,6 @@ namespace netgen
// which leads to an "order of destruction"-problem,
// thus we use the C-variant:
if (getenv ("NGPROFILE"))
{
char filename[100];

View File

@ -1,4 +1,4 @@
noinst_HEADERS = geom2dmesh.hpp spline2d.hpp splinegeometry.hpp spline.hpp geometry2d.hpp vsgeom2d.hpp
noinst_HEADERS = geom2dmesh.hpp splinegeometry.hpp spline.hpp geometry2d.hpp vsgeom2d.hpp
AM_CPPFLAGS = -I$(top_srcdir)/libsrc/include $(TCL_INCLUDES)

View File

@ -976,9 +976,9 @@ void SplineGeometry<D> :: LoadData ( ifstream & infile )
template<int D>
void SplineGeometry<D> :: PartitionBoundary (double h, Mesh & mesh2d)
void SplineGeometry2d :: PartitionBoundary (double h, Mesh & mesh2d)
{
enum { D = 2 };
Box<D> bbox;
GetBoundingBox (bbox);
double dist = Dist (bbox.PMin(), bbox.PMax());

View File

@ -38,6 +38,7 @@ namespace netgen
template < int D >
class SplineGeometry
{
protected:
Array < GeomPoint<D> > geompoints;
Array < SplineSeg<D>* > splines;
double elto0;
@ -48,7 +49,6 @@ namespace netgen
Array<bool> tensormeshing;
Array<int> layer;
private:
void AppendSegment(SplineSeg<D> * spline, const int leftdomain, const int rightdomain,
const int bc,
const double reffac, const bool hprefleft, const bool hprefright,
@ -65,7 +65,6 @@ namespace netgen
void LoadDataNew ( ifstream & infile );
void LoadDataV2 ( ifstream & infile );
void PartitionBoundary (double h, Mesh & mesh2d);
void GetRawData (Array<double> & raw_data) const;
@ -152,6 +151,8 @@ namespace netgen
virtual int GenerateMesh (Mesh*& mesh, MeshingParameters & mparam,
int perfstepsstart, int perfstepsend);
void PartitionBoundary (double h, Mesh & mesh2d);
virtual Refinement & GetRefinement () const;
};

View File

@ -27,9 +27,9 @@
#include <typeinfo>
#ifdef PARALLEL
#undef SEEK_SET
#undef SEEK_CUR
#undef SEEK_END
// #undef SEEK_SET
// #undef SEEK_CUR
// #undef SEEK_END
#include <mpi.h>
#endif

View File

@ -137,7 +137,7 @@ namespace netgen
for(int elind = 1; elind <= ne; elind++)
{
// Extract the current volume element
const Element & el = mesh.VolumeElement(elind);
// const Element & el = mesh.VolumeElement(elind);
// Get the face numbers of the faces of the current volume element
// The values returned are given a sign depending on the orientation
@ -596,7 +596,6 @@ namespace netgen
void WriteOpenFOAM15xFormat (const Mesh & mesh, const string & casename)
{
int i,j;
bool error = false;
char casefiles[256];
@ -627,8 +626,8 @@ namespace netgen
return;
}
if((mesh.SurfaceElement(nse/2).GetType() != TRIG)
&& (mesh.SurfaceElement(nse/2).GetType() != QUAD)
if(( (mesh.SurfaceElement(nse/2).GetType() != TRIG)
&& (mesh.SurfaceElement(nse/2).GetType() != QUAD) )
|| (mesh.VolumeElement(ne/2).GetType() == TET10)
|| (mesh.VolumeElement(ne/2).GetType() == PRISM12))
{

View File

@ -14,7 +14,9 @@ clusters.hpp hprefinement.hpp improve3.hpp meshtype.hpp \
METASOURCES = AUTO
noinst_LTLIBRARIES = libmesh.la
lib_LTLIBRARIES = libmesh.la
libmesh_la_SOURCES = adfront2.cpp adfront3.cpp bisect.cpp boundarylayer.cpp \
clusters.cpp curvedelems.cpp delaunay.cpp delaunay2d.cpp \
geomsearch.cpp global.cpp hprefinement.cpp improve2.cpp \

View File

@ -5,14 +5,17 @@
#ifdef METIS
namespace metis { extern "C" {
namespace metis
{
extern "C" {
#include <metis.h>
} }
}
}
using namespace metis;
#endif
using namespace metis;
namespace netgen
{
@ -311,8 +314,6 @@ namespace netgen
// distribute the mesh to the slave processors
// call it only for the master !
void Mesh :: Distribute ()
@ -326,7 +327,12 @@ namespace netgen
// partition mesh
#ifdef METIS
ParallelMetis ();
#else
for (ElementIndex ei = 0; ei < GetNE(); ei++)
(*this)[ei].SetPartition(ntasks * ei/GetNE() + 1);
#endif
#ifdef SCALASCA
#pragma pomp inst end(metis)
@ -349,7 +355,7 @@ namespace netgen
#ifdef METIS
void Mesh :: ParallelMetis ( )
{
int timer = NgProfiler::CreateTimer ("Mesh::Partition");
@ -475,12 +481,12 @@ namespace netgen
delete [] epart;
delete [] npart;
}
#endif
void Mesh :: PartHybridMesh () // Array<int> & neloc )
{
#ifdef METIS
int ne = GetNE();
int nn = GetNP();
@ -571,11 +577,15 @@ namespace netgen
delete [] xadj;
delete [] part;
delete [] adjacency;
#else
cout << "parthybridmesh not available" << endl;
#endif
}
void Mesh :: PartDualHybridMesh ( ) // Array<int> & neloc )
{
#ifdef METIS
int ne = GetNE();
// int nn = GetNP();
@ -677,6 +687,10 @@ namespace netgen
delete [] xadj;
delete [] part;
delete [] adjacency;
#else
cout << "partdualmesh not available" << endl;
#endif
}
@ -685,6 +699,7 @@ namespace netgen
void Mesh :: PartDualHybridMesh2D ( )
{
#ifdef METIS
int ne = GetNSE();
int nv = GetNV();
@ -758,6 +773,10 @@ namespace netgen
for (SurfaceElementIndex sei = 0; sei < ne; sei++)
(*this) [sei].SetPartition (part[sei]+1);
#else
cout << "partdualmesh not available" << endl;
#endif
}

View File

@ -7,6 +7,44 @@
namespace netgen
{
template <class T>
void QuickSortRec (FlatArray<T> & data,
int left, int right)
{
int i = left;
int j = right;
T midval = data[(left+right)/2];
do
{
while (data[i] < midval) i++;
while (midval < data[j]) j--;
if (i <= j)
{
Swap (data[i], data[j]);
i++; j--;
}
}
while (i <= j);
if (left < j) QuickSortRec (data, left, j);
if (i < right) QuickSortRec (data, i, right);
}
template <class T>
void QuickSort (FlatArray<T> & data)
{
if (data.Size() > 1)
QuickSortRec (data, 0, data.Size()-1);
}
MeshTopology :: MeshTopology (const Mesh & amesh)
: mesh(amesh)
{
@ -85,8 +123,7 @@ namespace netgen
for (ElementIndex ei = 0; ei < ne; ei++)
{
const Element & el = mesh[ei];
int nelv = el.GetNV();
for (int j = 0; j < nelv; j++)
for (int j = 0; j < el.GetNV(); j++)
cnt[el[j]]++;
}
@ -94,8 +131,7 @@ namespace netgen
for (ElementIndex ei = 0; ei < ne; ei++)
{
const Element & el = mesh[ei];
int nelv = el.GetNV();
for (int j = 0; j < nelv; j++)
for (int j = 0; j < el.GetNV(); j++)
vert2element->AddSave (el[j], ei+1);
}
@ -103,8 +139,7 @@ namespace netgen
for (SurfaceElementIndex sei = 0; sei < nse; sei++)
{
const Element2d & el = mesh[sei];
int nelv = el.GetNV();
for (int j = 0; j < nelv; j++)
for (int j = 0; j < el.GetNV(); j++)
cnt[el[j]]++;
}
@ -112,8 +147,7 @@ namespace netgen
for (SurfaceElementIndex sei = 0; sei < nse; sei++)
{
const Element2d & el = mesh[sei];
int nelv = el.GetNV();
for (int j = 0; j < nelv; j++)
for (int j = 0; j < el.GetNV(); j++)
vert2surfelement->AddSave (el[j], sei+1);
}
@ -134,7 +168,6 @@ namespace netgen
}
if (buildedges)
{
static int timer1 = NgProfiler::CreateTimer ("topology::buildedges");
@ -197,6 +230,8 @@ namespace netgen
Array<int,PointIndex::BASE> edgenr(nv);
Array<int,PointIndex::BASE> edgeflag(nv);
Array<int> vertex2;
edgeflag = 0;
ned = edge2vert.Size();
@ -204,6 +239,8 @@ namespace netgen
for (int i = PointIndex::BASE; i < nv+PointIndex::BASE; i++)
{
vertex2.SetSize (0);
for (int j = 0; j < vert2edge[i].Size(); j++)
{
int ednr = vert2edge[i][j];
@ -211,15 +248,95 @@ namespace netgen
edgeflag[i2] = i;
edgenr[i2] = ednr;
}
for (int j = 0; j < vert2vertcoarse[i].Size(); j++) // fix by Markus
{
int v2 = vert2vertcoarse[i][j];
if (edgeflag[v2] < i)
{
ned++;
edgenr[v2] = ned;
*testout << "do we really need that ??????????????" << endl;
// ned++;
// edgenr[v2] = ned;
edgeflag[v2] = i;
missing.Append (INDEX_3(i,v2,ned));
vertex2.Append (v2);
// missing.Append (INDEX_3(i,v2,ned));
}
}
for (int j = 0; j < (*vert2element)[i].Size(); j++)
{
int elnr = (*vert2element)[i][j];
const Element & el = mesh.VolumeElement (elnr);
int neledges = GetNEdges (el.GetType());
const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType());
for (int k = 0; k < neledges; k++)
{
INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]);
edge.Sort();
if (edge.I1() != i) continue;
if (edgeflag[edge.I2()] < i)
{
vertex2.Append (edge.I2());
edgeflag[edge.I2()] = i;
}
}
}
for (int j = 0; j < (*vert2surfelement)[i].Size(); j++)
{
int elnr = (*vert2surfelement)[i][j];
const Element2d & el = mesh.SurfaceElement (elnr);
int neledges = GetNEdges (el.GetType());
const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType());
for (int k = 0; k < neledges; k++)
{
INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]);
edge.Sort();
if (edge.I1() != i) continue;
if (edgeflag[edge.I2()] < i)
{
vertex2.Append (edge.I2());
edgeflag[edge.I2()] = i;
}
}
}
for (int j = 0; j < (*vert2segment)[i].Size(); j++)
{
int elnr = (*vert2segment)[i][j];
const Segment & el = mesh.LineSegment (elnr);
INDEX_2 edge(el[0], el[1]);
edge.Sort();
if (edge.I1() != i) continue;
if (edgeflag[edge.I2()] < i)
{
vertex2.Append (edge.I2());
edgeflag[edge.I2()] = i;
}
}
QuickSort (vertex2);
for (int j = 0; j < vertex2.Size(); j++)
edgenr[vertex2[j]] = ++ned;
for (int j = 0; j < vert2vertcoarse[i].Size(); j++) // fix by Markus
{
int v2 = vert2vertcoarse[i][j];
if (edgeflag[v2] < i)
{
// ned++;
// edgenr[v2] = ned;
// edgeflag[v2] = i;
missing.Append (INDEX_3(i,v2,edgenr[v2]));
}
}
@ -241,13 +358,6 @@ namespace netgen
if (edge.I1() != i)
continue;
if (edgeflag[edge.I2()] < i)
{
ned++;
edgenr[edge.I2()] = ned;
edgeflag[edge.I2()] = i;
}
int edgenum = edgenr[edge.I2()];
if (edgedir) edgenum *= -1;
edges.Elem(elnr)[k] = edgenum;
@ -269,15 +379,7 @@ namespace netgen
int edgedir = (edge.I1() > edge.I2());
if (edgedir) swap (edge.I1(), edge.I2());
if (edge.I1() != i)
continue;
if (edgeflag[edge.I2()] < i)
{
ned++;
edgenr[edge.I2()] = ned;
edgeflag[edge.I2()] = i;
}
if (edge.I1() != i) continue;
int edgenum = edgenr[edge.I2()];
if (edgedir) edgenum *= -1;
@ -295,23 +397,16 @@ namespace netgen
int edgedir = (edge.I1() > edge.I2());
if (edgedir) swap (edge.I1(), edge.I2());
if (edge.I1() != i)
continue;
if (edge.I1() != i) continue;
if (edgeflag[edge.I2()] < i)
{
ned++;
edgenr[edge.I2()] = ned;
edgeflag[edge.I2()] = i;
}
int edgenum = edgenr[edge.I2()];
if (edgedir) edgenum *= -1;
segedges.Elem(elnr) = edgenum;
}
}
edge2vert.SetSize (ned);
for (int i = 1; i <= ne; i++)
{
@ -333,6 +428,11 @@ namespace netgen
edge2vert.Elem(edgenum)[1] = edge.I2();
}
}
*testout << "edge 2 vert:" << endl;
for (int i = 0; i < edge2vert.Size(); i++)
*testout << edge2vert[i][0] << " " << edge2vert[i][1] << endl;
for (int i = 1; i <= nse; i++)
{
const Element2d & el = mesh.SurfaceElement (i);

View File

@ -752,7 +752,7 @@ namespace netgen
GLXContextID xid = glXGetContextIDEXT (ctx);
displname = XDisplayName (0);
/*
cout << "Init Parallel GL" << endl;
cout << "DisplayName = " << displname << endl;
cout << "current display = " << dpy << endl;
@ -762,6 +762,7 @@ namespace netgen
cout << "contextid = " << xid << endl;
cout << "isdirect = " << glXIsDirect ( dpy, ctx ) << endl;
cout << "extensionstring = " << glXQueryExtensionsString( dpy, 0 ) << endl;
*/
Array<MPI_Request> request(ntasks);
MPI_Status status;

View File

@ -1020,6 +1020,7 @@ namespace netgen
#ifdef PARALLEL
/*
if ( el.IsGhost() )
{
if ( faceindex == selface )
@ -1028,6 +1029,7 @@ namespace netgen
glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat_coll_transp);
}
else
*/
{
if ( faceindex == selface )
glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, matcolsel);
@ -1894,9 +1896,20 @@ namespace netgen
for (ElementIndex ei = 0; ei < mesh->GetNE(); ei++)
{
if (vispar.drawtetsdomain > 0 &&
vispar.drawtetsdomain != (*mesh)[ei].GetIndex())
continue;
if (vispar.drawtetsdomain > 0)
{
int tetid =
#ifdef PARALLEL
vispar.drawmetispartition ?
(*mesh)[ei].GetPartition()
:
#endif
(*mesh)[ei].GetIndex();
if (vispar.drawtetsdomain != tetid)
continue;
}
const Element & el = (*mesh)[ei];

View File

@ -12,8 +12,6 @@ netgen_LDADD = $(top_builddir)/libsrc/visualization/libvisual.a \
$(top_builddir)/libsrc/geom2d/libgeom2dvis.la \
$(top_builddir)/libsrc/geom2d/libgeom2d.la \
$(top_builddir)/libsrc/interface/libinterface.la \
$(top_builddir)/libsrc/occ/liboccvis.la \
$(top_builddir)/libsrc/occ/libocc.la \
$(top_builddir)/libsrc/stlgeom/libstlvis.la \
$(top_builddir)/libsrc/stlgeom/libstl.la \
$(top_builddir)/libsrc/meshing/libmesh.la \
@ -22,6 +20,10 @@ netgen_LDADD = $(top_builddir)/libsrc/visualization/libvisual.a \
$(top_builddir)/libsrc/general/libgen.la \
$(OCCLIBS) -L$(TK_BIN_DIR)/Togl1.7 $(TOGLLIBDIR) -lTogl1.7 $(LIBGLU) $(TK_LIB_SPEC) $(TCL_LIB_SPEC) $(MPI_LIBS) $(FFMPEG_LIBS) $(JPEGLIB_LIBS) $(PKG_LIBS)
# $(top_builddir)/libsrc/occ/liboccvis.la
# $(top_builddir)/libsrc/occ/libocc.la
# add for static linkage of ngsolve:
# /opt/netgen/lib/libngsolve.a /opt/netgen/lib/libngcomp.a /opt/netgen/lib/libngcomp.a /opt/netgen/lib/libngfemng.a /opt/netgen/lib/libngmg.a /opt/netgen/lib/libngla.a /opt/netgen/lib/libngbla.a /opt/netgen/lib/libngstd.a -L/opt/intel/mkl/10.2.1.017/lib/em64t /opt/intel/mkl/10.2.1.017/lib/em64t/libmkl_solver_lp64.a -lmkl_intel_lp64 -lmkl_gnu_thread -lmkl_core
#

View File

@ -129,8 +129,8 @@ int main(int argc, char ** argv)
#ifdef PARALLEL
cout << "Including MPI " << endl;
cout << "Using " << netgen::ntasks << " processor"
cout << "Running MPI - parallel using "
<< netgen::ntasks << " processor"
<< ((netgen::ntasks > 1) ? "s " : " ") << endl;
#endif
}

View File

@ -1,5 +1,7 @@
#ifdef PARALLEL
#include "dlfcn.h"
#ifdef OCCGEOMETRY
#include <occgeom.hpp>
#endif
@ -29,7 +31,6 @@
#endif
#include <geometry2d.hpp>
#include <stlgeom.hpp>
@ -48,11 +49,8 @@
#include "parallelfunc.hpp"
// extern "C" void NGS_ParallelRun (const string & message);
void (*NGS_ParallelRun) (const string & message) = NULL;
@ -71,49 +69,17 @@ void Parallel_Exit();
namespace netgen {
extern AutoPtr<Mesh> mesh;
// geometry: either CSG, or, if an other is non-null,
// then the other
extern AutoPtr<CSGeometry> geometry ;
extern STLGeometry * stlgeometry;
extern AutoPtr<SplineGeometry2d> geometry2d ;
#ifdef OCCGEOMETRY
extern OCCGeometry * occgeometry;
#endif
#ifdef ACIS
extern ACISGeometry * acisgeometry;
#endif
}
using namespace netgen;
using netgen::RegisterUserFormats;
#ifdef PARALLEL
void Ng_Exit ()
{
#ifdef NGSOLVE
Parallel_Exit();
#endif
delete stlgeometry;
stlgeometry = NULL;
#ifdef OCCGEOMETRY
delete occgeometry;
occgeometry = 0;
#endif
geometry.Reset (NULL);
geometry2d.Reset (NULL);
// delete testout;
return;
// Parallel_Exit();
}
#endif
@ -152,15 +118,31 @@ void ParallelRun()
#pragma pomp inst end (message)
#endif
#ifdef NGSOLVE
if ( message.compare(0, 3, "ngs") == 0 )
{
PrintMessage ( 2, "Starting NgSolve routine ", message ) ;
NGS_ParallelRun (message);
PrintMessage ( 1, "Starting NgSolve routine ", message ) ;
if (NGS_ParallelRun == NULL)
{
void * handle = dlopen ("libngsolve.so", RTLD_NOW | RTLD_GLOBAL);
if (!handle)
{
cerr << "cannot load shared library libngsolve.so" << endl;
exit(1);
}
NGS_ParallelRun = (void (*) (const string & message)) dlsym (handle, "NGS_ParallelRun");
if (!NGS_ParallelRun)
{
cerr << "cannot bind function NGS_ParallelRun" << endl;
exit(1);
}
}
(*NGS_ParallelRun) (message);
}
else
#endif
if ( message == "mesh" )
{
@ -274,14 +256,16 @@ void ParallelRun()
Window win;
int wx, wy;
unsigned int ww, wh, bw, depth;
cout << "got drawable: " << curDrawable << endl;
// cout << "got drawable: " << curDrawable << endl;
XGetGeometry(display, curDrawable, &win,
&wx, &wy, &ww, &wh,
&bw, &depth);
/*
cout << "P" << id << ": window-props: x = " << wx << ", y = " << wy
<< ", w = " << ww << ", h = " << wh << ", depth = " << depth << endl;
*/
#define VISUAL
#ifdef VISUAL
@ -348,7 +332,7 @@ void ParallelRun()
// context = glXCreateContext( display, visinfo, 0, /* curContext, */ False );
context = glXCreateContext( display, visinfo, glXImportContextEXT ( display, contextid ), False );
cout << "context = " << context << endl;
// cout << "context = " << context << endl;
glXMakeCurrent (display, curDrawable, context);
@ -387,7 +371,7 @@ void ParallelRun()
if (redraw_cmd == "filledlist")
{
glXMakeCurrent (display, curDrawable, context);
vsmesh.BuildFilledList();
vsmesh.BuildFilledList (0);
glXMakeCurrent (display, None, NULL);
}

View File

@ -10,7 +10,6 @@ namespace netgen
extern AutoPtr<Mesh> mesh;
using namespace netgen;
// int NgPar_Glob2Loc_SurfEl ( int globnum )
@ -32,7 +31,9 @@ namespace netgen
// {
// return mesh->GetParallelTopology().Glob2Loc_Vert(globnum+1) -1;
// }
}
using namespace netgen;
int NgPar_GetLoc2Glob_VolEl ( int locnum )
{
@ -152,6 +153,6 @@ namespace netgen
return mesh -> GetParallelTopology().IsGhostEdge ( edgenum+1);
}
}
#endif