summaryrefslogtreecommitdiff
path: root/benchmarks/CUDA/DG/src
diff options
context:
space:
mode:
Diffstat (limited to 'benchmarks/CUDA/DG/src')
-rw-r--r--benchmarks/CUDA/DG/src/BuildMaps2d.c231
-rw-r--r--benchmarks/CUDA/DG/src/BuildMaps3d.c279
-rw-r--r--benchmarks/CUDA/DG/src/FacePair2d.c154
-rw-r--r--benchmarks/CUDA/DG/src/FacePair3d.c201
-rw-r--r--benchmarks/CUDA/DG/src/InitCPU3d.c152
-rw-r--r--benchmarks/CUDA/DG/src/LoadBalance3d.c243
-rw-r--r--benchmarks/CUDA/DG/src/MaxwellsDriver3d.c110
-rw-r--r--benchmarks/CUDA/DG/src/MaxwellsKernel3d.cu551
-rw-r--r--benchmarks/CUDA/DG/src/MaxwellsMPI3d.c71
-rw-r--r--benchmarks/CUDA/DG/src/MaxwellsRHS3d.c220
-rw-r--r--benchmarks/CUDA/DG/src/MaxwellsRun3d.c60
-rw-r--r--benchmarks/CUDA/DG/src/Mesh3d.c195
-rw-r--r--benchmarks/CUDA/DG/src/ParallelPairs.c124
-rw-r--r--benchmarks/CUDA/DG/src/StartUp3d.c123
-rw-r--r--benchmarks/CUDA/DG/src/Utils.c155
-rw-r--r--benchmarks/CUDA/DG/src/buildsource.m106
-rw-r--r--benchmarks/CUDA/DG/src/gpuDG.vpj173
-rw-r--r--benchmarks/CUDA/DG/src/tictoc.cu23
18 files changed, 3171 insertions, 0 deletions
diff --git a/benchmarks/CUDA/DG/src/BuildMaps2d.c b/benchmarks/CUDA/DG/src/BuildMaps2d.c
new file mode 100644
index 0000000..7d0d57c
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/BuildMaps2d.c
@@ -0,0 +1,231 @@
+#include "mpi.h"
+#include <fem.h>
+
+void BuildMaps2d(Mesh *mesh){
+
+ int nprocs = mesh->nprocs;
+ int procid = mesh->procid;
+
+ int K = mesh->K;
+ int Nfaces = mesh->Nfaces;
+
+ mesh->vmapM = BuildIntVector(p_Nfp*p_Nfaces*K);
+ mesh->vmapP = BuildIntVector(p_Nfp*p_Nfaces*K);
+
+ int m;
+ int k1,f1,p1,n1,id1, k2,f2,p2,n2,id2;
+
+ double x1, y1, x2, y2, d12;
+
+ double *nxk = BuildVector(Nfaces);
+ double *nyk = BuildVector(Nfaces);
+ double *sJk = BuildVector(Nfaces);
+
+ /* first build local */
+ for(k1=0;k1<K;++k1){
+
+ /* get some information about the face geometries */
+ Normals2d(mesh, k1, nxk, nyk, sJk);
+
+ for(f1=0;f1<Nfaces;++f1){
+
+ /* volume -> face nodes */
+ for(n1=0;n1<p_Nfp;++n1){
+ id1 = n1+f1*p_Nfp+k1*p_Nfp*p_Nfaces;
+ mesh->vmapM[id1] = mesh->Fmask[f1][n1] + k1*p_Np;
+ }
+
+
+ /* find neighbor */
+ k2 = mesh->EToE[k1][f1];
+ f2 = mesh->EToF[k1][f1];
+ p2 = mesh->EToP[k1][f1];
+
+ if(k1==k2 || procid!=p2 ){
+ for(n1=0;n1<p_Nfp;++n1){
+ id1 = n1+f1*p_Nfp+k1*p_Nfp*p_Nfaces;
+ mesh->vmapP[id1] = k1*p_Np + mesh->Fmask[f1][n1];
+ }
+ }else{
+ /* treat as boundary for the moment */
+
+ for(n1=0;n1<p_Nfp;++n1){
+ id1 = n1+f1*p_Nfp+k1*p_Nfp*p_Nfaces;
+
+ x1 = mesh->x[k1][mesh->Fmask[f1][n1]];
+ y1 = mesh->y[k1][mesh->Fmask[f1][n1]];
+ for(n2=0;n2<p_Nfp;++n2){
+
+ id2 = n2+f2*p_Nfp+k2*p_Nfp*p_Nfaces;
+
+ x2 = mesh->x[k2][mesh->Fmask[f2][n2]];
+ y2 = mesh->y[k2][mesh->Fmask[f2][n2]];
+
+ /* find normalized distance between these nodes */
+ /* [ use sJk as a measure of edge length (ignore factor of 2) ] */
+ d12 = ((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2))/(sJk[f1]*sJk[f1]);
+ if(d12<NODETOL){
+ mesh->vmapP[id1] = k2*p_Np + mesh->Fmask[f2][n2];
+ }
+ }
+ }
+ }
+ }
+ }
+
+#if 0
+ int n;
+ for(n=0;n<p_Nfp*p_Nfaces*mesh->K;++n){
+ x1 = mesh->x[0][mesh->vmapM[n]];
+ y1 = mesh->y[0][mesh->vmapM[n]];
+ x2 = mesh->x[0][mesh->vmapP[n]];
+ y2 = mesh->y[0][mesh->vmapP[n]];
+ d12 = ((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
+ printf("n:%d %d -> %d d=%lg\n", n, mesh->vmapM[n], mesh->vmapP[n], d12);
+ }
+#endif
+
+ /* now build parallel maps */
+ double **xsend = (double**) calloc(nprocs, sizeof(double*));
+ double **ysend = (double**) calloc(nprocs, sizeof(double*));
+ double **xrecv = (double**) calloc(nprocs, sizeof(double*));
+ double **yrecv = (double**) calloc(nprocs, sizeof(double*));
+
+ int **Esend = (int**) calloc(nprocs, sizeof(int*));
+ int **Fsend = (int**) calloc(nprocs, sizeof(int*));
+ int **Erecv = (int**) calloc(nprocs, sizeof(int*));
+ int **Frecv = (int**) calloc(nprocs, sizeof(int*));
+
+ for(p2=0;p2<nprocs;++p2){
+ if(mesh->Npar[p2]){
+ xsend[p2] = BuildVector(mesh->Npar[p2]*p_Nfp);
+ ysend[p2] = BuildVector(mesh->Npar[p2]*p_Nfp);
+ Esend[p2] = BuildIntVector(mesh->Npar[p2]*p_Nfp);
+ Fsend[p2] = BuildIntVector(mesh->Npar[p2]*p_Nfp);
+
+ xrecv[p2] = BuildVector(mesh->Npar[p2]*p_Nfp);
+ yrecv[p2] = BuildVector(mesh->Npar[p2]*p_Nfp);
+ Erecv[p2] = BuildIntVector(mesh->Npar[p2]*p_Nfp);
+ Frecv[p2] = BuildIntVector(mesh->Npar[p2]*p_Nfp);
+ }
+ }
+
+ int *skP = BuildIntVector(nprocs);
+
+ /* send coordinates in local order */
+ int cnt = 0;
+ for(k1=0;k1<K;++k1){
+ for(f1=0;f1<p_Nfaces;++f1){
+ p2 = mesh->EToP[k1][f1];
+ if(p2!=procid){
+ for(n1=0;n1<p_Nfp;++n1){
+ xsend[p2][skP[p2]] = mesh->x[k1][mesh->Fmask[f1][n1]];
+ ysend[p2][skP[p2]] = mesh->y[k1][mesh->Fmask[f1][n1]];
+ Esend[p2][skP[p2]] = mesh->EToE[k1][f1];
+ Fsend[p2][skP[p2]] = mesh->EToF[k1][f1];
+ ++(skP[p2]);
+ }
+ }
+ }
+ }
+
+ MPI_Request *xsendrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *ysendrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *xrecvrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *yrecvrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *Esendrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *Fsendrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *Erecvrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *Frecvrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+
+ MPI_Status *status = (MPI_Status*) calloc(nprocs, sizeof(MPI_Status));
+
+ cnt = 0;
+ for(p2=0;p2<nprocs;++p2){
+ if(p2!=procid && mesh->Npar[p2]!=0){
+ int Nout = mesh->Npar[p2]*p_Nfp;
+
+ MPI_Isend(xsend[p2], Nout, MPI_DOUBLE, p2, 666+p2, MPI_COMM_WORLD, xsendrequests+cnt);
+ MPI_Isend(ysend[p2], Nout, MPI_DOUBLE, p2, 1666+p2, MPI_COMM_WORLD, ysendrequests+cnt);
+ MPI_Isend(Esend[p2], Nout, MPI_INT, p2, 2666+p2, MPI_COMM_WORLD, Esendrequests+cnt);
+ MPI_Isend(Fsend[p2], Nout, MPI_INT, p2, 3666+p2, MPI_COMM_WORLD, Fsendrequests+cnt);
+
+ MPI_Irecv(xrecv[p2], Nout, MPI_DOUBLE, p2, 666+procid, MPI_COMM_WORLD, xrecvrequests+cnt);
+ MPI_Irecv(yrecv[p2], Nout, MPI_DOUBLE, p2, 1666+procid, MPI_COMM_WORLD, yrecvrequests+cnt);
+ MPI_Irecv(Erecv[p2], Nout, MPI_INT, p2, 2666+procid, MPI_COMM_WORLD, Erecvrequests+cnt);
+ MPI_Irecv(Frecv[p2], Nout, MPI_INT, p2, 3666+procid, MPI_COMM_WORLD, Frecvrequests+cnt);
+ ++cnt;
+ }
+ }
+
+ MPI_Waitall(cnt, xsendrequests, status);
+ MPI_Waitall(cnt, ysendrequests, status);
+ MPI_Waitall(cnt, Esendrequests, status);
+ MPI_Waitall(cnt, Fsendrequests, status);
+ MPI_Waitall(cnt, xrecvrequests, status);
+ MPI_Waitall(cnt, yrecvrequests, status);
+ MPI_Waitall(cnt, Erecvrequests, status);
+ MPI_Waitall(cnt, Frecvrequests, status);
+
+ /* add up the total number of outgoing/ingoing nodes */
+ mesh->parNtotalout = 0;
+ for(p2=0;p2<nprocs;++p2)
+ mesh->parNtotalout += skP[p2]*p_Nfields;
+
+ mesh->parmapOUT = BuildIntVector(mesh->parNtotalout);
+
+ /* now match up local nodes with the requested (recv'ed nodes) */
+ int idout = -1;
+ int sk = 0;
+ for(p2=0;p2<nprocs;++p2){
+ /* for each received face */
+ for(m=0;m<skP[p2];++m){
+ k1 = Erecv[p2][m];
+ f1 = Frecv[p2][m];
+ x2 = xrecv[p2][m];
+ y2 = yrecv[p2][m];
+
+ Normals2d(mesh, k1, nxk, nyk, sJk);
+
+ for(n1=0;n1<p_Nfp;++n1){
+
+ x1 = mesh->x[k1][mesh->Fmask[f1][n1]];
+ y1 = mesh->y[k1][mesh->Fmask[f1][n1]];
+
+ d12 = ((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2))/(sJk[f1]*sJk[f1]);
+
+ if(d12<NODETOL){
+ int fld;
+ for(fld=0;fld<p_Nfields;++fld){
+#ifdef CUDA
+ mesh->parmapOUT[sk++] = k1*BSIZE*p_Nfields+mesh->Fmask[f1][n1] + BSIZE*fld;
+#else
+ mesh->parmapOUT[sk++] = p_Nfields*(k1*p_Np+mesh->Fmask[f1][n1]) + fld;
+#endif
+ }
+ }
+ }
+ }
+ }
+
+ /* create incoming node map */
+ int parcnt = -1;
+ for(p2=0;p2<nprocs;++p2){
+ for(k1=0;k1<K;++k1){
+ for(f1=0;f1<p_Nfaces;++f1){
+ if(mesh->EToP[k1][f1]==p2 && p2!=procid){
+ for(n1=0;n1<p_Nfp;++n1){
+ id1 = n1+f1*p_Nfp+k1*p_Nfp*p_Nfaces;
+ mesh->vmapP[id1] = parcnt;
+ --parcnt;
+ }
+ }
+ }
+ }
+ }
+
+ /* buffers for communication */
+ mesh->f_outQ = (float*) calloc(mesh->parNtotalout+1, sizeof(float));
+ mesh->f_inQ = (float*) calloc(mesh->parNtotalout+1, sizeof(float));
+
+}
diff --git a/benchmarks/CUDA/DG/src/BuildMaps3d.c b/benchmarks/CUDA/DG/src/BuildMaps3d.c
new file mode 100644
index 0000000..89f7a44
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/BuildMaps3d.c
@@ -0,0 +1,279 @@
+#include "mpi.h"
+#include "fem.h"
+
+void BuildMaps3d(Mesh *mesh){
+
+ printf("Hello %d\n", 1002);
+
+ int nprocs = mesh->nprocs;
+ int procid = mesh->procid;
+
+ int K = mesh->K;
+ int Nfaces = mesh->Nfaces;
+
+ mesh->vmapM = BuildIntVector(p_Nfp*p_Nfaces*K);
+ mesh->vmapP = BuildIntVector(p_Nfp*p_Nfaces*K);
+
+ int m;
+ int k1,f1,p1,n1,id1, k2,f2,p2,n2,id2;
+
+ double x1, y1,z1, x2, y2, z2, d12;
+
+ double *nxk = BuildVector(Nfaces);
+ double *nyk = BuildVector(Nfaces);
+ double *nzk = BuildVector(Nfaces);
+ double *sJk = BuildVector(Nfaces);
+
+ printf("Hello %d\n", 1001);
+
+ /* first build local */
+ for(k1=0;k1<K;++k1){
+
+ /* get some information about the face geometries */
+ Normals3d(mesh, k1, nxk, nyk, nzk, sJk);
+
+ for(f1=0;f1<Nfaces;++f1){
+
+ /* volume -> face nodes */
+ for(n1=0;n1<p_Nfp;++n1){
+ id1 = n1+f1*p_Nfp+k1*p_Nfp*p_Nfaces;
+ mesh->vmapM[id1] = mesh->Fmask[f1][n1] + k1*p_Np;
+ }
+
+
+ /* find neighbor */
+ k2 = mesh->EToE[k1][f1];
+ f2 = mesh->EToF[k1][f1];
+ p2 = mesh->EToP[k1][f1];
+
+ if(k1==k2 || procid!=p2 ){
+ for(n1=0;n1<p_Nfp;++n1){
+ id1 = n1+f1*p_Nfp+k1*p_Nfp*p_Nfaces;
+ mesh->vmapP[id1] = k1*p_Np + mesh->Fmask[f1][n1];
+ }
+ }else{
+ /* treat as boundary for the moment */
+
+ for(n1=0;n1<p_Nfp;++n1){
+ id1 = n1+f1*p_Nfp+k1*p_Nfp*p_Nfaces;
+
+ x1 = mesh->x[k1][mesh->Fmask[f1][n1]];
+ y1 = mesh->y[k1][mesh->Fmask[f1][n1]];
+ z1 = mesh->z[k1][mesh->Fmask[f1][n1]];
+
+ for(n2=0;n2<p_Nfp;++n2){
+
+ id2 = n2+f2*p_Nfp+k2*p_Nfp*p_Nfaces;
+
+ x2 = mesh->x[k2][mesh->Fmask[f2][n2]];
+ y2 = mesh->y[k2][mesh->Fmask[f2][n2]];
+ z2 = mesh->z[k2][mesh->Fmask[f2][n2]];
+
+ /* find normalized distance between these nodes */
+ /* [ use sJk as a measure of edge length (ignore factor of 2) ] */
+ d12 = ((x1-x2)*(x1-x2) +
+ (y1-y2)*(y1-y2) +
+ (z1-z2)*(z1-z2)); /* /(sJk[f1]*sJk[f1]); */
+ if(d12<NODETOL){
+ mesh->vmapP[id1] = k2*p_Np + mesh->Fmask[f2][n2];
+ break;
+ }
+ }
+ if(n2==p_Nfp){
+ printf("LOST NODE !!!\n");
+ }
+ }
+ }
+ }
+ }
+
+#if 0
+ int n;
+ for(k1=0;k1<mesh->K;++k1){
+ double drdx, dsdx, dtdx;
+ double drdy, dsdy, dtdy;
+ double drdz, dsdz, dtdz, J;
+
+ Normals3d(mesh, k1, nxk, nyk, nzk, sJk);
+
+ GeometricFactors3d(mesh, k1,
+ &drdx, &dsdx, &dtdx,
+ &drdy, &dsdy, &dtdy,
+ &drdz, &dsdz, &dtdz, &J);
+
+ for(f1=0;f1<mesh->Nfaces;++f1){
+ for(m=0;m<p_Nfp;++m){
+ n = k1*p_Nfp*p_Nfaces+f1*p_Nfp+m;
+ x1 = mesh->x[0][mesh->vmapM[n]];
+ y1 = mesh->y[0][mesh->vmapM[n]];
+ z1 = mesh->z[0][mesh->vmapM[n]];
+ x2 = mesh->x[0][mesh->vmapP[n]];
+ y2 = mesh->y[0][mesh->vmapP[n]];
+ z2 = mesh->z[0][mesh->vmapP[n]];
+ d12 = ((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2) );
+ printf("n:%d %d -> %d d=%lg sJ=%lg J=%lg (%d,%d,%d,%d)\n",
+ n, mesh->vmapM[n], mesh->vmapP[n], d12, sJk[f1], J,
+ mesh->EToV[k1][0],mesh->EToV[k1][1],mesh->EToV[k1][2],mesh->EToV[k1][3]);
+ }
+ }
+ }
+#endif
+
+ /* now build parallel maps */
+ double **xsend = (double**) calloc(nprocs, sizeof(double*));
+ double **ysend = (double**) calloc(nprocs, sizeof(double*));
+ double **zsend = (double**) calloc(nprocs, sizeof(double*));
+ double **xrecv = (double**) calloc(nprocs, sizeof(double*));
+ double **yrecv = (double**) calloc(nprocs, sizeof(double*));
+ double **zrecv = (double**) calloc(nprocs, sizeof(double*));
+
+ int **Esend = (int**) calloc(nprocs, sizeof(int*));
+ int **Fsend = (int**) calloc(nprocs, sizeof(int*));
+ int **Erecv = (int**) calloc(nprocs, sizeof(int*));
+ int **Frecv = (int**) calloc(nprocs, sizeof(int*));
+
+ for(p2=0;p2<nprocs;++p2){
+ if(mesh->Npar[p2]){
+ xsend[p2] = BuildVector(mesh->Npar[p2]*p_Nfp);
+ ysend[p2] = BuildVector(mesh->Npar[p2]*p_Nfp);
+ zsend[p2] = BuildVector(mesh->Npar[p2]*p_Nfp);
+ Esend[p2] = BuildIntVector(mesh->Npar[p2]*p_Nfp);
+ Fsend[p2] = BuildIntVector(mesh->Npar[p2]*p_Nfp);
+
+ xrecv[p2] = BuildVector(mesh->Npar[p2]*p_Nfp);
+ yrecv[p2] = BuildVector(mesh->Npar[p2]*p_Nfp);
+ zrecv[p2] = BuildVector(mesh->Npar[p2]*p_Nfp);
+ Erecv[p2] = BuildIntVector(mesh->Npar[p2]*p_Nfp);
+ Frecv[p2] = BuildIntVector(mesh->Npar[p2]*p_Nfp);
+ }
+ }
+
+ int *skP = BuildIntVector(nprocs);
+
+ /* send coordinates in local order */
+ int cnt = 0;
+ for(k1=0;k1<K;++k1){
+ for(f1=0;f1<p_Nfaces;++f1){
+ p2 = mesh->EToP[k1][f1];
+ if(p2!=procid){
+ for(n1=0;n1<p_Nfp;++n1){
+ xsend[p2][skP[p2]] = mesh->x[k1][mesh->Fmask[f1][n1]];
+ ysend[p2][skP[p2]] = mesh->y[k1][mesh->Fmask[f1][n1]];
+ zsend[p2][skP[p2]] = mesh->z[k1][mesh->Fmask[f1][n1]];
+ Esend[p2][skP[p2]] = mesh->EToE[k1][f1];
+ Fsend[p2][skP[p2]] = mesh->EToF[k1][f1];
+ ++(skP[p2]);
+ }
+ }
+ }
+ }
+
+ MPI_Request *xsendrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *ysendrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *zsendrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *xrecvrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *yrecvrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *zrecvrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *Esendrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *Fsendrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *Erecvrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *Frecvrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+
+ MPI_Status *status = (MPI_Status*) calloc(nprocs, sizeof(MPI_Status));
+
+ cnt = 0;
+ for(p2=0;p2<nprocs;++p2){
+ if(p2!=procid && mesh->Npar[p2]!=0){
+ int Nout = mesh->Npar[p2]*p_Nfp;
+
+ MPI_Isend(xsend[p2], Nout, MPI_DOUBLE, p2, 666+p2, MPI_COMM_WORLD, xsendrequests+cnt);
+ MPI_Isend(ysend[p2], Nout, MPI_DOUBLE, p2, 1666+p2, MPI_COMM_WORLD, ysendrequests+cnt);
+ MPI_Isend(zsend[p2], Nout, MPI_DOUBLE, p2, 4666+p2, MPI_COMM_WORLD, zsendrequests+cnt);
+ MPI_Isend(Esend[p2], Nout, MPI_INT, p2, 2666+p2, MPI_COMM_WORLD, Esendrequests+cnt);
+ MPI_Isend(Fsend[p2], Nout, MPI_INT, p2, 3666+p2, MPI_COMM_WORLD, Fsendrequests+cnt);
+
+ MPI_Irecv(xrecv[p2], Nout, MPI_DOUBLE, p2, 666+procid, MPI_COMM_WORLD, xrecvrequests+cnt);
+ MPI_Irecv(yrecv[p2], Nout, MPI_DOUBLE, p2, 1666+procid, MPI_COMM_WORLD, yrecvrequests+cnt);
+ MPI_Irecv(zrecv[p2], Nout, MPI_DOUBLE, p2, 4666+procid, MPI_COMM_WORLD, zrecvrequests+cnt);
+ MPI_Irecv(Erecv[p2], Nout, MPI_INT, p2, 2666+procid, MPI_COMM_WORLD, Erecvrequests+cnt);
+ MPI_Irecv(Frecv[p2], Nout, MPI_INT, p2, 3666+procid, MPI_COMM_WORLD, Frecvrequests+cnt);
+ ++cnt;
+ }
+ }
+
+ MPI_Waitall(cnt, xsendrequests, status);
+ MPI_Waitall(cnt, ysendrequests, status);
+ MPI_Waitall(cnt, zsendrequests, status);
+ MPI_Waitall(cnt, Esendrequests, status);
+ MPI_Waitall(cnt, Fsendrequests, status);
+
+ MPI_Waitall(cnt, xrecvrequests, status);
+ MPI_Waitall(cnt, yrecvrequests, status);
+ MPI_Waitall(cnt, zrecvrequests, status);
+ MPI_Waitall(cnt, Erecvrequests, status);
+ MPI_Waitall(cnt, Frecvrequests, status);
+
+ /* add up the total number of outgoing/ingoing nodes */
+ mesh->parNtotalout = 0;
+ for(p2=0;p2<nprocs;++p2)
+ mesh->parNtotalout += skP[p2]*p_Nfields;
+
+ mesh->parmapOUT = BuildIntVector(mesh->parNtotalout);
+
+ /* now match up local nodes with the requested (recv'ed nodes) */
+ int idout = -1;
+ int sk = 0;
+ for(p2=0;p2<nprocs;++p2){
+ /* for each received face */
+ for(m=0;m<skP[p2];++m){
+ k1 = Erecv[p2][m];
+ f1 = Frecv[p2][m];
+ x2 = xrecv[p2][m];
+ y2 = yrecv[p2][m];
+ z2 = zrecv[p2][m];
+
+ Normals3d(mesh, k1, nxk, nyk, nzk, sJk);
+
+ for(n1=0;n1<p_Nfp;++n1){
+
+ x1 = mesh->x[k1][mesh->Fmask[f1][n1]];
+ y1 = mesh->y[k1][mesh->Fmask[f1][n1]];
+ z1 = mesh->z[k1][mesh->Fmask[f1][n1]];
+
+ d12 = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2)); /* /(sJk[f1]*sJk[f1]); */
+
+ if(d12<NODETOL){
+ int fld;
+ for(fld=0;fld<p_Nfields;++fld){
+#ifdef CUDA
+ mesh->parmapOUT[sk++] = k1*BSIZE*p_Nfields+mesh->Fmask[f1][n1] + BSIZE*fld;
+#else
+ mesh->parmapOUT[sk++] = p_Nfields*(k1*p_Np+mesh->Fmask[f1][n1]) + fld;
+#endif
+ }
+ }
+ }
+ }
+ }
+
+ /* create incoming node map */
+ int parcnt = -1;
+ for(p2=0;p2<nprocs;++p2){
+ for(k1=0;k1<K;++k1){
+ for(f1=0;f1<p_Nfaces;++f1){
+ if(mesh->EToP[k1][f1]==p2 && p2!=procid){
+ for(n1=0;n1<p_Nfp;++n1){
+ id1 = n1+f1*p_Nfp+k1*p_Nfp*p_Nfaces;
+ mesh->vmapP[id1] = parcnt;
+ --parcnt;
+ }
+ }
+ }
+ }
+ }
+
+ /* buffers for communication */
+ mesh->f_outQ = (float*) calloc(mesh->parNtotalout+1, sizeof(float));
+ mesh->f_inQ = (float*) calloc(mesh->parNtotalout+1, sizeof(float));
+
+}
diff --git a/benchmarks/CUDA/DG/src/FacePair2d.c b/benchmarks/CUDA/DG/src/FacePair2d.c
new file mode 100644
index 0000000..92e47f5
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/FacePair2d.c
@@ -0,0 +1,154 @@
+#include "mpi.h"
+#include "fem.h"
+
+typedef struct foob {
+ int p1, k1, f1, p2, k2, f2, va, vb, g;
+}face;
+
+int compare_pairs(const void *obj1, const void *obj2){
+
+ face *e1 = (face*) obj1;
+ face *e2 = (face*) obj2;
+
+ int a1 = e1->va, b1 = e1->vb;
+ int a2 = e2->va, b2 = e2->vb;
+
+ int va1, vb1, va2, vb2;
+
+ va1 = min(a1, b1);
+ vb1 = max(a1, b1);
+
+ va2 = min(a2, b2);
+ vb2 = max(a2, b2);
+
+ if(vb1<vb2)
+ return -1;
+ else if(vb1>vb2)
+ return 1;
+ else if(va1<va2)
+ return -1;
+ else if(va1>va2)
+ return 1;
+
+ return 0;
+}
+
+int pairprocget(const void *obj1){
+ face *e1 = (face*) obj1;
+ return (e1->p1);
+}
+
+
+int pairnumget(const void *obj1){
+ face *e1 = (face*) obj1;
+ return (e1->g);
+}
+
+void pairnumset(const void *obj1, int g){
+ face *e1 = (face*) obj1;
+ e1->g = g;
+}
+
+void pairmarry(const void *obj1, const void *obj2){
+
+ face *e1 = (face*) obj1;
+ face *e2 = (face*) obj2;
+ e1->p2 = e2->p1; e1->k2 = e2->k1; e1->f2 = e2->f1;
+ e2->p2 = e1->p1; e2->k2 = e1->k1; e2->f2 = e1->f1;
+}
+
+void FacePair2d(Mesh *mesh, int *maxNv){
+
+ int procid = mesh->procid;
+ int nprocs = mesh->nprocs;
+
+ int Klocal = mesh->K;
+ int Nfaces = mesh->Nfaces;
+ int Nverts = mesh->Nverts;
+
+ int **EToV = mesh->EToV;
+
+ const int vnum[3][2] = { {0,1}, {1,2}, {2,0} };
+
+ int n, k, e, sk, v;
+
+ face *myfaces = (face*) calloc(Klocal*Nfaces, sizeof(face));
+
+ /* find maximum local vertex number */
+ int localmaxgnum = 0;
+ for(k=0;k<Klocal;++k)
+ for(v=0;v<Nverts;++v)
+ localmaxgnum = max(localmaxgnum, EToV[k][v]);
+ ++localmaxgnum;
+
+ int maxgnum;
+ MPI_Allreduce(&localmaxgnum, &maxgnum, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
+ *maxNv = maxgnum;
+
+ sk = 0;
+ for(k=0;k<Klocal;++k){
+ for(e=0;e<Nfaces;++e){
+ int n1 = EToV[k][vnum[e][0]];
+ int n2 = EToV[k][vnum[e][1]];
+
+ myfaces[sk].p1 = procid; myfaces[sk].k1 = k; myfaces[sk].f1 = e;
+ myfaces[sk].p2 = procid; myfaces[sk].k2 = k; myfaces[sk].f2 = e;
+
+ myfaces[sk].va = max(n1,n2);
+ myfaces[sk].vb = min(n1,n2);
+ myfaces[sk].g = max(n1,n2); /* marker for sorting into bins */
+ ++sk;
+ }
+ }
+
+ ParallelPairs(myfaces, Klocal*Nfaces, sizeof(face),
+ pairnumget, pairnumset, pairprocget,
+ pairmarry, compare_pairs);
+
+ mesh->Npar = BuildIntVector(nprocs);
+
+ mesh->EToE = BuildIntMatrix(Klocal, Nfaces);
+ mesh->EToF = BuildIntMatrix(Klocal, Nfaces);
+ mesh->EToP = BuildIntMatrix(Klocal, Nfaces);
+
+ int id, k1, k2, f1, f2, p1, p2;
+ sk = 0;
+
+ for(n=0;n<Klocal*Nfaces;++n){
+
+ k1 = myfaces[n].k1; f1 = myfaces[n].f1; p1 = myfaces[n].p1;
+ k2 = myfaces[n].k2; f2 = myfaces[n].f2; p2 = myfaces[n].p2;
+
+ if(p1!=procid)
+ fprintf(stderr, "WARNING WRONG proc\n");
+
+ mesh->EToE[k1][f1] = k2;
+ mesh->EToF[k1][f1] = f2;
+ mesh->EToP[k1][f1] = p2;
+
+ if(p1!=p2){
+ /* increment number of links */
+ ++mesh->Npar[p2];
+ }
+ }
+
+ mesh->parK = (int**) calloc(nprocs, sizeof(int*));
+ mesh->parF = (int**) calloc(nprocs, sizeof(int*));
+ for(p2=0;p2<nprocs;++p2){
+ mesh->parK[p2] = BuildIntVector(mesh->Npar[p2]);
+ mesh->parF[p2] = BuildIntVector(mesh->Npar[p2]);
+ mesh->Npar[p2] = 0;
+ for(n=0;n<Klocal*Nfaces;++n){
+ if(myfaces[n].p2==p2 && p2!=procid){
+ int k1 = myfaces[n].k1, f1 = myfaces[n].f1;
+ int k2 = myfaces[n].k2, f2 = myfaces[n].f2;
+ mesh->parK[p2][mesh->Npar[p2] ] = k1;
+ mesh->parF[p2][mesh->Npar[p2]++] = f1;
+ }
+ }
+ }
+
+ free(myfaces);
+}
+
+
diff --git a/benchmarks/CUDA/DG/src/FacePair3d.c b/benchmarks/CUDA/DG/src/FacePair3d.c
new file mode 100644
index 0000000..41d1fad
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/FacePair3d.c
@@ -0,0 +1,201 @@
+#include "mpi.h"
+#include "fem.h"
+
+typedef struct foob {
+ int p1, k1, f1, p2, k2, f2;
+ int va, vb, vc, g;
+}face3d;
+
+int compare_pairs3d(const void *obj1, const void *obj2){
+
+ face3d *e1 = (face3d*) obj1;
+ face3d *e2 = (face3d*) obj2;
+
+ int a1 = e1->va, b1 = e1->vb, c1 = e1->vc;
+ int a2 = e2->va, b2 = e2->vb, c2 = e2->vc;
+
+ int va1, vb1, vc1, va2, vb2, vc2;
+
+ va1 = min(a1, min(b1, c1));
+ vc1 = max(a1, max(b1, c1));
+
+ if(va1!=a1 && vc1!=a1) vb1=a1;
+ if(va1!=b1 && vc1!=b1) vb1=b1;
+ if(va1!=c1 && vc1!=c1) vb1=c1;
+
+ va2 = min(a2, min(b2, c2));
+ vc2 = max(a2, max(b2, c2));
+
+ if(va2!=a2 && vc2!=a2) vb2=a2;
+ if(va2!=b2 && vc2!=b2) vb2=b2;
+ if(va2!=c2 && vc2!=c2) vb2=c2;
+
+ if(vc1<vc2)
+ return -1;
+ else if(vc1>vc2)
+ return 1;
+ else if(vb1<vb2)
+ return -1;
+ else if(vb1>vb2)
+ return 1;
+ else if(va1<va2)
+ return -1;
+ else if(va1>va2)
+ return 1;
+
+ return 0;
+
+}
+
+
+int pairprocget3d(const void *obj1){
+ face3d *e1 = (face3d*) obj1;
+ return (e1->p1);
+}
+
+
+int pairnumget3d(const void *obj1){
+ face3d *e1 = (face3d*) obj1;
+ return (e1->g);
+}
+
+void pairnumset3d(const void *obj1, int g){
+ face3d *e1 = (face3d*) obj1;
+ e1->g = g;
+}
+
+void pairmarry3d(const void *obj1, const void *obj2){
+
+ face3d *e1 = (face3d*) obj1;
+ face3d *e2 = (face3d*) obj2;
+ e1->p2 = e2->p1; e1->k2 = e2->k1; e1->f2 = e2->f1;
+ e2->p2 = e1->p1; e2->k2 = e1->k1; e2->f2 = e1->f1;
+}
+
+void FacePair3d(Mesh *mesh, int *maxNv){
+
+ int procid = mesh->procid;
+ int nprocs = mesh->nprocs;
+
+ int Klocal = mesh->K;
+ int Nfaces = mesh->Nfaces;
+ int Nverts = mesh->Nverts;
+
+ int **EToV = mesh->EToV;
+
+ const int vnum[4][3] = { {0,1,2}, {0,1,3}, {1,2,3}, {0,2,3} };
+
+ int n, k, e, sk, v;
+
+ face3d *myfaces = (face3d*) calloc(Klocal*Nfaces, sizeof(face3d));
+
+ /* find maximum local vertex number */
+ int localmaxgnum = 0;
+ for(k=0;k<Klocal;++k)
+ for(v=0;v<Nverts;++v)
+ localmaxgnum = max(localmaxgnum, EToV[k][v]);
+ ++localmaxgnum;
+
+ int maxgnum;
+ MPI_Allreduce(&localmaxgnum, &maxgnum, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
+ *maxNv = maxgnum;
+
+ sk = 0;
+ for(k=0;k<Klocal;++k){
+ for(e=0;e<Nfaces;++e){
+ int a1 = EToV[k][vnum[e][0]];
+ int b1 = EToV[k][vnum[e][1]];
+ int c1 = EToV[k][vnum[e][2]];
+
+ myfaces[sk].p1 = procid; myfaces[sk].k1 = k; myfaces[sk].f1 = e;
+ myfaces[sk].p2 = procid; myfaces[sk].k2 = k; myfaces[sk].f2 = e;
+
+ int va1, vb1, vc1, va2, vb2, vc2;
+
+ va1 = min(a1, min(b1, c1));
+ vc1 = max(a1, max(b1, c1));
+
+ if(va1!=a1 && vc1!=a1) vb1=a1;
+ if(va1!=b1 && vc1!=b1) vb1=b1;
+ if(va1!=c1 && vc1!=c1) vb1=c1;
+
+ myfaces[sk].va = va1;
+ myfaces[sk].vb = vb1;
+ myfaces[sk].vc = vc1;
+ myfaces[sk].g = max(va1,max(vb1,vc1)); /* marker for sorting into bins */
+ ++sk;
+ }
+ }
+
+ ParallelPairs(myfaces, Klocal*Nfaces, sizeof(face3d),
+ pairnumget3d, pairnumset3d, pairprocget3d,
+ pairmarry3d, compare_pairs3d);
+
+ mesh->Npar = BuildIntVector(nprocs);
+
+ mesh->EToE = BuildIntMatrix(Klocal, Nfaces);
+ mesh->EToF = BuildIntMatrix(Klocal, Nfaces);
+ mesh->EToP = BuildIntMatrix(Klocal, Nfaces);
+
+ int id, k1, k2, f1, f2, p1, p2;
+ sk = 0;
+
+ for(n=0;n<Klocal*Nfaces;++n){
+
+ k1 = myfaces[n].k1; f1 = myfaces[n].f1; p1 = myfaces[n].p1;
+ k2 = myfaces[n].k2; f2 = myfaces[n].f2; p2 = myfaces[n].p2;
+
+ if(p1!=procid){
+ fprintf(stderr, "WARNING WRONG proc\n");
+ exit(-1);
+ }
+
+ mesh->EToE[k1][f1] = k2;
+ mesh->EToF[k1][f1] = f2;
+ mesh->EToP[k1][f1] = p2;
+
+ if(p1!=p2){
+ /* increment number of links */
+ ++mesh->Npar[p2];
+ }
+ }
+
+#if 0
+ char fname[BUFSIZ];
+ sprintf(fname, "proc%d.dat", mesh->procid);
+ FILE *fp = fopen(fname, "w");
+ for(k1=0;k1<mesh->K;++k1){
+ for(f1=0;f1<mesh->Nfaces;++f1){
+ fprintf(fp, "p: (%d,%d,%d)->(%d,%d,%d)\n",
+ k1,f1,mesh->procid,
+ mesh->EToE[k1][f1],
+ mesh->EToF[k1][f1],
+ mesh->EToP[k1][f1]);
+
+ }
+ }
+ fclose(fp);
+#endif
+
+ mesh->parK = (int**) calloc(nprocs, sizeof(int*));
+ mesh->parF = (int**) calloc(nprocs, sizeof(int*));
+ for(p2=0;p2<nprocs;++p2){
+ mesh->parK[p2] = BuildIntVector(mesh->Npar[p2]);
+ mesh->parF[p2] = BuildIntVector(mesh->Npar[p2]);
+ mesh->Npar[p2] = 0;
+ for(n=0;n<Klocal*Nfaces;++n){
+ if(myfaces[n].p2==p2 && p2!=procid){
+ int k1 = myfaces[n].k1, f1 = myfaces[n].f1;
+ int k2 = myfaces[n].k2, f2 = myfaces[n].f2;
+ mesh->parK[p2][mesh->Npar[p2] ] = k1;
+ mesh->parF[p2][mesh->Npar[p2]++] = f1;
+ }
+ }
+ printf("proc: %d sends %d to proc: %d\n",
+ mesh->procid, mesh->Npar[p2], p2);
+ }
+
+ free(myfaces);
+}
+
+
diff --git a/benchmarks/CUDA/DG/src/InitCPU3d.c b/benchmarks/CUDA/DG/src/InitCPU3d.c
new file mode 100644
index 0000000..2fe907f
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/InitCPU3d.c
@@ -0,0 +1,152 @@
+#include "mpi.h"
+#include "fem.h"
+
+double InitCPU3d(Mesh *mesh, int Nfields){
+
+ printf("Np = %d, BSIZE = %d\n", p_Np, BSIZE);
+
+ /* Q */
+ int sz = mesh->K*(p_Np)*Nfields*sizeof(float); /* TW BLOCK */
+
+ mesh->f_Q = (float*) calloc(mesh->K*p_Np*Nfields, sizeof(float));
+ mesh->f_rhsQ = (float*) calloc(mesh->K*p_Np*Nfields, sizeof(float));
+ mesh->f_resQ = (float*) calloc(mesh->K*p_Np*Nfields, sizeof(float));
+
+ /* float LIFT */
+ sz = p_Np*(p_Nfp)*(p_Nfaces)*sizeof(float);
+ mesh->f_LIFT = (float*) malloc(sz);
+ int sk = 0, n, m, f, k;
+
+ for(n=0;n<p_Np;++n){
+ for(m=0;m<p_Nfp*p_Nfaces;++m){
+ mesh->f_LIFT[sk++] = mesh->LIFT[n][m];
+ }
+ }
+
+ /* float Dr & Ds */
+ sz = p_Np*p_Np*sizeof(float);
+ mesh->f_Dr = (float*) malloc(sz);
+ mesh->f_Ds = (float*) malloc(sz);
+ mesh->f_Dt = (float*) malloc(sz);
+
+ sk = 0;
+ for(n=0;n<p_Np;++n){
+ for(m=0;m<p_Np;++m){
+ mesh->f_Dr[sk] = mesh->Dr[n][m];
+ mesh->f_Ds[sk] = mesh->Ds[n][m];
+ mesh->f_Dt[sk] = mesh->Dt[n][m];
+ ++sk;
+ }
+ }
+
+ /* vgeo */
+ double drdx, dsdx, dtdx;
+ double drdy, dsdy, dtdy;
+ double drdz, dsdz, dtdz, J;
+ mesh->vgeo = (float*) calloc(12*mesh->K, sizeof(float));
+
+ for(k=0;k<mesh->K;++k){
+ GeometricFactors3d(mesh, k,
+ &drdx, &dsdx, &dtdx,
+ &drdy, &dsdy, &dtdy,
+ &drdz, &dsdz, &dtdz, &J);
+
+ mesh->vgeo[k*12+0] = drdx; mesh->vgeo[k*12+1] = drdy; mesh->vgeo[k*12+2] = drdz;
+ mesh->vgeo[k*12+4] = dsdx; mesh->vgeo[k*12+5] = dsdy; mesh->vgeo[k*12+6] = dsdz;
+ mesh->vgeo[k*12+8] = dtdx; mesh->vgeo[k*12+9] = dtdy; mesh->vgeo[k*12+10] = dtdz;
+ }
+
+ /* surfinfo (vmapM, vmapP, Fscale, Bscale, nx, ny, nz, 0) */
+ sz = mesh->K*p_Nfp*p_Nfaces*7*sizeof(float);
+
+ mesh->surfinfo = (float*) malloc(sz);
+
+ /* local-local info */
+ sk = 0;
+ int skP = -1;
+ double *nxk = BuildVector(mesh->Nfaces);
+ double *nyk = BuildVector(mesh->Nfaces);
+ double *nzk = BuildVector(mesh->Nfaces);
+ double *sJk = BuildVector(mesh->Nfaces);
+
+ double dt = 1e6;
+
+ sk = 0;
+ for(k=0;k<mesh->K;++k){
+
+ GeometricFactors3d(mesh, k,
+ &drdx, &dsdx, &dtdx,
+ &drdy, &dsdy, &dtdy,
+ &drdz, &dsdz, &dtdz, &J);
+
+ Normals3d(mesh, k, nxk, nyk, nzk, sJk);
+
+ for(f=0;f<mesh->Nfaces;++f){
+
+ dt = min(dt, J/sJk[f]);
+
+ for(m=0;m<p_Nfp;++m){
+ int id = m + f*p_Nfp + p_Nfp*p_Nfaces*k;
+ int idM = mesh->vmapM[id];
+ int idP = mesh->vmapP[id];
+ int nM = idM%p_Np;
+ int nP = idP%p_Np;
+ int kM = (idM-nM)/p_Np;
+ int kP = (idP-nP)/p_Np;
+ idM = Nfields*(nM + p_Np*kM);
+ idP = Nfields*(nP + p_Np*kP);
+
+ /* stub resolve some other way */
+ if(mesh->vmapP[id]<0){
+ idP = mesh->vmapP[id]; /* -ve numbers */
+ }
+
+ mesh->surfinfo[sk++] = idM;
+ mesh->surfinfo[sk++] = idP;
+ mesh->surfinfo[sk++] = sJk[f]/(2.*J);
+ mesh->surfinfo[sk++] = (idM==idP)?-1.:1.;
+ mesh->surfinfo[sk++] = nxk[f];
+ mesh->surfinfo[sk++] = nyk[f];
+ mesh->surfinfo[sk++] = nzk[f];
+ }
+ }
+ }
+}
+
+void cpu_set_data3d(Mesh *mesh, double *Hx, double *Hy, double *Hz,
+ double *Ex, double *Ey, double *Ez){
+
+ const int K = mesh->K;
+ int k, n, sk=0;
+
+ for(k=0;k<K;++k){
+ for(n=0;n<p_Np;++n){
+ mesh->f_Q[sk++] = Hx[n+k*p_Np];
+ mesh->f_Q[sk++] = Hy[n+k*p_Np];
+ mesh->f_Q[sk++] = Hz[n+k*p_Np];
+ mesh->f_Q[sk++] = Ex[n+k*p_Np];
+ mesh->f_Q[sk++] = Ey[n+k*p_Np];
+ mesh->f_Q[sk++] = Ez[n+k*p_Np];
+ }
+ }
+
+}
+
+
+void cpu_get_data3d(Mesh *mesh, double *Hx, double *Hy, double *Hz,
+ double *Ex, double *Ey, double *Ez){
+ const int K = mesh->K;
+ int k, n, sk=0;
+
+ for(k=0;k<K;++k){
+ for(n=0;n<p_Np;++n){
+ Hx[n+k*p_Np] = mesh->f_Q[sk++];
+ Hy[n+k*p_Np] = mesh->f_Q[sk++];
+ Hz[n+k*p_Np] = mesh->f_Q[sk++];
+ Ex[n+k*p_Np] = mesh->f_Q[sk++];
+ Ey[n+k*p_Np] = mesh->f_Q[sk++];
+ Ez[n+k*p_Np] = mesh->f_Q[sk++];
+ }
+ }
+
+}
diff --git a/benchmarks/CUDA/DG/src/LoadBalance3d.c b/benchmarks/CUDA/DG/src/LoadBalance3d.c
new file mode 100644
index 0000000..586cd8c
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/LoadBalance3d.c
@@ -0,0 +1,243 @@
+#include "mpi.h"
+#include <parmetisbin.h>
+
+#include "fem.h"
+
+void LoadBalance3d(Mesh *mesh){
+
+ int n,p,k,v,f;
+
+ int nprocs = mesh->nprocs;
+ int procid = mesh->procid;
+ int **EToV = mesh->EToV;
+ double **VX = mesh->GX;
+ double **VY = mesh->GY;
+ double **VZ = mesh->GZ;
+
+ if(!procid) printf("Root: Entering LoadBalance\n");
+
+ int Nverts = mesh->Nverts;
+
+ int *Kprocs = BuildIntVector(nprocs);
+
+ /* local number of elements */
+ int Klocal = mesh->K;
+
+ /* find number of elements on all processors */
+ MPI_Allgather(&Klocal, 1, MPI_INT, Kprocs, 1, MPI_INT, MPI_COMM_WORLD);
+
+ /* element distribution -- cumulative element count on processes */
+ idxtype *elmdist = idxmalloc(nprocs+1, "elmdist");
+
+ elmdist[0] = 0;
+ for(p=0;p<nprocs;++p)
+ elmdist[p+1] = elmdist[p] + Kprocs[p];
+
+ /* list of element starts */
+ idxtype *eptr = idxmalloc(Klocal+1, "eptr");
+
+ eptr[0] = 0;
+ for(k=0;k<Klocal;++k)
+ eptr[k+1] = eptr[k] + Nverts;
+
+ /* local element to vertex */
+ idxtype *eind = idxmalloc(Nverts*Klocal, "eind");
+
+ for(k=0;k<Klocal;++k)
+ for(n=0;n<Nverts;++n)
+ eind[k*Nverts+n] = EToV[k][n];
+
+ /* weight per element */
+ idxtype *elmwgt = idxmalloc(Klocal, "elmwgt");
+
+ for(k=0;k<Klocal;++k)
+ elmwgt[k] = 1.;
+
+ /* weight flag */
+ int wgtflag = 0;
+
+ /* number flag (1=fortran, 0=c) */
+ int numflag = 0;
+
+ /* ncon = 1 */
+ int ncon = 1;
+
+ /* nodes on element face */
+ int ncommonnodes = 3;
+
+ /* number of partitions */
+ int nparts = nprocs;
+
+ /* tpwgts */
+ float *tpwgts = (float*) calloc(Klocal, sizeof(float));
+
+ for(k=0;k<Klocal;++k)
+ tpwgts[k] = 1./(float)nprocs;
+
+ float ubvec[MAXNCON];
+
+ for (n=0; n<ncon; ++n)
+ ubvec[n] = UNBALANCE_FRACTION;
+
+ int options[10];
+
+ options[0] = 1;
+ options[PMV3_OPTION_DBGLVL] = 7;
+ options[PMV3_OPTION_SEED] = 0;
+
+ int edgecut;
+
+ idxtype *part = idxmalloc(Klocal, "part");
+
+ MPI_Comm comm;
+ MPI_Comm_dup(MPI_COMM_WORLD, &comm);
+
+ ParMETIS_V3_PartMeshKway
+ (elmdist,
+ eptr,
+ eind,
+ elmwgt,
+ &wgtflag,
+ &numflag,
+ &ncon,
+ &ncommonnodes,
+ &nparts,
+ tpwgts,
+ ubvec,
+ options,
+ &edgecut,
+ part,
+ &comm);
+
+ int **outlist = (int**) calloc(nprocs, sizeof(int*));
+ double **xoutlist = (double**) calloc(nprocs, sizeof(double*));
+ double **youtlist = (double**) calloc(nprocs, sizeof(double*));
+ double **zoutlist = (double**) calloc(nprocs, sizeof(double*));
+
+ int *outK = (int*) calloc(nprocs, sizeof(int));
+
+ int *inK = (int*) calloc(nprocs, sizeof(int));
+
+ MPI_Request *inrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *outrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *xinrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *xoutrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *yinrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *youtrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *zinrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ MPI_Request *zoutrequests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+
+ for(k=0;k<Klocal;++k)
+ ++outK[part[k]];
+
+ /* get count of incoming elements from each process */
+ MPI_Alltoall(outK, 1, MPI_INT,
+ inK, 1, MPI_INT,
+ MPI_COMM_WORLD);
+
+ /* count totals on each process */
+ int * newKprocs = BuildIntVector(nprocs);
+ MPI_Allreduce(outK, newKprocs, nprocs, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+
+ int totalinK = 0;
+ for(p=0;p<nprocs;++p){
+ totalinK += inK[p];
+ }
+
+ int **newEToV = BuildIntMatrix(totalinK, Nverts);
+ double **newVX = BuildMatrix(totalinK, Nverts);
+ double **newVY = BuildMatrix(totalinK, Nverts);
+ double **newVZ = BuildMatrix(totalinK, Nverts);
+
+ int cnt = 0;
+ for(p=0;p<nprocs;++p){
+ MPI_Irecv(newEToV[cnt], Nverts*inK[p], MPI_INT, p, 666+p, MPI_COMM_WORLD,
+ inrequests+p);
+ MPI_Irecv(newVX[cnt], Nverts*inK[p], MPI_DOUBLE, p, 1666+p, MPI_COMM_WORLD,
+ xinrequests+p);
+ MPI_Irecv(newVY[cnt], Nverts*inK[p], MPI_DOUBLE, p, 2666+p, MPI_COMM_WORLD,
+ yinrequests+p);
+ MPI_Irecv(newVZ[cnt], Nverts*inK[p], MPI_DOUBLE, p, 3666+p, MPI_COMM_WORLD,
+ zinrequests+p);
+ cnt = cnt + inK[p];
+ }
+
+ for(p=0;p<nprocs;++p){
+ int cnt = 0;
+ outlist[p] = BuildIntVector(Nverts*outK[p]);
+ xoutlist[p] = BuildVector(Nverts*outK[p]);
+ youtlist[p] = BuildVector(Nverts*outK[p]);
+ zoutlist[p] = BuildVector(Nverts*outK[p]);
+
+ for(k=0;k<Klocal;++k)
+ if(part[k]==p){
+ for(v=0;v<Nverts;++v){
+ outlist[p][cnt] = EToV[k][v];
+ xoutlist[p][cnt] = VX[k][v];
+ youtlist[p][cnt] = VY[k][v];
+ zoutlist[p][cnt] = VZ[k][v];
+ ++cnt;
+ }
+ }
+
+ MPI_Isend(outlist[p], Nverts*outK[p], MPI_INT, p, 666+procid, MPI_COMM_WORLD,
+ outrequests+p);
+ MPI_Isend(xoutlist[p], Nverts*outK[p], MPI_DOUBLE, p, 1666+procid, MPI_COMM_WORLD,
+ xoutrequests+p);
+ MPI_Isend(youtlist[p], Nverts*outK[p], MPI_DOUBLE, p, 2666+procid, MPI_COMM_WORLD,
+ youtrequests+p);
+ MPI_Isend(zoutlist[p], Nverts*outK[p], MPI_DOUBLE, p, 3666+procid, MPI_COMM_WORLD,
+ zoutrequests+p);
+ }
+
+ MPI_Status *instatus = (MPI_Status*) calloc(nprocs, sizeof(MPI_Status));
+ MPI_Status *outstatus = (MPI_Status*) calloc(nprocs, sizeof(MPI_Status));
+
+ MPI_Waitall(nprocs, inrequests, instatus);
+ MPI_Waitall(nprocs, xinrequests, instatus);
+ MPI_Waitall(nprocs, yinrequests, instatus);
+ MPI_Waitall(nprocs, zinrequests, instatus);
+
+ MPI_Waitall(nprocs, outrequests, outstatus);
+ MPI_Waitall(nprocs, xoutrequests, outstatus);
+ MPI_Waitall(nprocs, youtrequests, outstatus);
+ MPI_Waitall(nprocs, zoutrequests, outstatus);
+
+ if(mesh->GX!=NULL){
+ DestroyMatrix(mesh->GX);
+ DestroyMatrix(mesh->GY);
+ DestroyMatrix(mesh->GZ);
+ DestroyIntMatrix(mesh->EToV);
+ }
+
+ mesh->GX = newVX;
+ mesh->GY = newVY;
+ mesh->GZ = newVZ;
+ mesh->EToV = newEToV;
+ mesh->K = totalinK;
+
+ for(p=0;p<nprocs;++p){
+ if(outlist[p]){
+ free(outlist[p]);
+ free(xoutlist[p]);
+ free(youtlist[p]);
+ free(zoutlist[p]);
+ }
+ }
+
+ free(outK);
+ free(inK);
+
+ free(inrequests);
+ free(outrequests);
+
+ free(xinrequests);
+ free(xoutrequests);
+ free(yinrequests);
+ free(youtrequests);
+ free(zinrequests);
+ free(zoutrequests);
+ free(instatus);
+ free(outstatus);
+
+}
diff --git a/benchmarks/CUDA/DG/src/MaxwellsDriver3d.c b/benchmarks/CUDA/DG/src/MaxwellsDriver3d.c
new file mode 100644
index 0000000..dd761c7
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/MaxwellsDriver3d.c
@@ -0,0 +1,110 @@
+#include "mpi.h"
+#include "fem.h"
+
+main(int argc, char **argv){
+
+ Mesh *mesh;
+ int procid, nprocs, maxNv;
+ int k,n, sk=0;
+ double minEz, maxEz, gminEz, gmaxEz;
+
+ /* initialize MPI */
+ MPI_Init(&argc, &argv);
+
+ /* assign gpu */
+ MPI_Comm_rank(MPI_COMM_WORLD, &procid);
+ MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+
+#ifdef CUDA
+ cudaSetDevice((procid+2)%4);
+#endif
+
+ /* (parallel) read part of fem mesh from file */
+ mesh = ReadMesh3d(argv[1]);
+
+ /* perform load balancing */
+ LoadBalance3d(mesh);
+
+ /* find element-element connectivity */
+ FacePair3d(mesh, &maxNv);
+
+ /* perform start up */
+ StartUp3d(mesh);
+
+ /* field storage (double) */
+ double *Hx = (double*) calloc(mesh->K*p_Np, sizeof(double));
+ double *Hy = (double*) calloc(mesh->K*p_Np, sizeof(double));
+ double *Hz = (double*) calloc(mesh->K*p_Np, sizeof(double));
+ double *Ex = (double*) calloc(mesh->K*p_Np, sizeof(double));
+ double *Ey = (double*) calloc(mesh->K*p_Np, sizeof(double));
+ double *Ez = (double*) calloc(mesh->K*p_Np, sizeof(double));
+
+ /* initial conditions */
+ for(k=0;k<mesh->K;++k){
+ for(n=0;n<p_Np;++n) {
+ Hx[sk] = 0;
+ Hy[sk] = 0;
+ Hz[sk] = 0;
+ Ex[sk] = 0;
+ Ey[sk] = 0;
+ Ez[sk] = cos(M_PI*mesh->x[k][n])*cos(M_PI*mesh->y[k][n])*cos(M_PI*mesh->z[k][n]);
+ ++sk;
+ }
+ }
+
+ double dt, gdt;
+
+#ifdef CUDA
+ /* initialize GPU info */
+ dt = InitGPU3d(mesh, p_Nfields);
+
+ /* load data onto GPU */
+ gpu_set_data3d(mesh->K, Hx, Hy, Hz, Ex, Ey, Ez);
+#else
+ /* initialize GPU info */
+ dt = InitCPU3d(mesh, p_Nfields);
+
+ /* load data onto CPU float storage */
+ cpu_set_data3d(mesh, Hx, Hy, Hz, Ex, Ey, Ez);
+#endif
+
+ MPI_Allreduce(&dt, &gdt, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+
+ dt = .5*gdt/((p_N+1)*(p_N+1));
+
+ if(mesh->procid==0)
+ printf("dt = %f\n", dt);
+
+ double FinalTime = .00050;
+
+ /* solve */
+ MaxwellsRun3d(mesh, FinalTime, dt);
+
+#ifdef CUDA
+ /* unload data from GPU */
+ gpu_get_data3d(mesh->K, Hx, Hy, Hz, Ex, Ey, Ez);
+#else
+ cpu_get_data3d(mesh, Hx, Hy, Hz, Ex, Ey, Ez);
+#endif
+
+ /* find maximum & minimum values for Ez */
+ minEz=Ez[0], maxEz=Ez[0];
+
+ for(n=0;n<mesh->K*p_Np;++n) {
+ minEz = (minEz>Ez[n])?Ez[n]:minEz;
+ maxEz = (maxEz<Ez[n])?Ez[n]:maxEz;
+ printf("%e\n", Ez[n]);
+ }
+
+ MPI_Reduce(&minEz, &gminEz, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD);
+ MPI_Reduce(&maxEz, &gmaxEz, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
+
+ if(procid==0)
+ printf("t=%f Ez in [ %f, %f ] \n", FinalTime, gminEz, gmaxEz );
+
+ /* nicely stop MPI */
+ MPI_Finalize();
+
+ /* end game */
+ exit(1);
+}
diff --git a/benchmarks/CUDA/DG/src/MaxwellsKernel3d.cu b/benchmarks/CUDA/DG/src/MaxwellsKernel3d.cu
new file mode 100644
index 0000000..bfa4396
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/MaxwellsKernel3d.cu
@@ -0,0 +1,551 @@
+/* -*- mode: C; c-basic-offset: 8; c-indent-level: 8; c-continued-statement-offset: 8; c-label-offset: -8; -*- */
+
+#include <stdio.h>
+#include <cuda.h>
+
+texture<float4, 1, cudaReadModeElementType> t_LIFT;
+texture<float4, 1, cudaReadModeElementType> t_DrDsDt;
+texture<float, 1, cudaReadModeElementType> t_Dr;
+texture<float, 1, cudaReadModeElementType> t_Ds;
+texture<float, 1, cudaReadModeElementType> t_Dt;
+texture<float, 1, cudaReadModeElementType> t_vgeo;
+texture<float4, 1, cudaReadModeElementType> t_vgeo4;
+texture<float, 1, cudaReadModeElementType> t_Q;
+texture<float, 1, cudaReadModeElementType> t_partQ;
+texture<float, 1, cudaReadModeElementType> t_surfinfo;
+
+static float *c_LIFT;
+static float *c_DrDsDt;
+static float *c_surfinfo;
+static float *c_vgeo;
+static float *c_Q;
+static float *c_partQ;
+static float *c_rhsQ;
+static float *c_resQ;
+static float *c_tmp;
+
+extern "C"
+{
+
+#include "fem.h"
+
+double InitGPU3d(Mesh *mesh, int Nfields){
+
+ /* Q */
+ int sz = mesh->K*(BSIZE)*p_Nfields*sizeof(float);
+
+ float *f_Q = (float*) calloc(mesh->K*BSIZE*p_Nfields, sizeof(float));
+ cudaMalloc ((void**) &c_Q, sz);
+ cudaMalloc ((void**) &c_rhsQ, sz);
+ cudaMalloc ((void**) &c_resQ, sz);
+ cudaMalloc ((void**) &c_tmp, sz);
+ cudaMemcpy( c_Q, f_Q, sz, cudaMemcpyHostToDevice);
+ cudaMemcpy( c_rhsQ, f_Q, sz, cudaMemcpyHostToDevice);
+ cudaMemcpy( c_resQ, f_Q, sz, cudaMemcpyHostToDevice);
+ cudaMemcpy( c_tmp, f_Q, sz, cudaMemcpyHostToDevice);
+
+ cudaBindTexture(0, t_Q, c_Q, sz);
+
+ sz = mesh->parNtotalout*sizeof(float);
+ cudaMalloc((void**) &c_partQ, sz);
+ cudaBindTexture(0, t_partQ, c_partQ, sz);
+
+ /* LIFT */
+ sz = p_Np*(p_Nfp)*p_Nfaces*sizeof(float);
+#if 0
+ float *f_LIFT = (float*) malloc(sz);
+ int skL = 0;
+ for(int m=0;m<p_Nfp*p_Nfaces;++m){
+ for(int n=0;n<p_Np;++n){
+ f_LIFT[skL++] = d_LIFT[n+p_Np*m];
+ }
+ }
+#else
+ float *f_LIFT = (float*) malloc(sz);
+ int skL = 0;
+ for(int m=0;m<p_Nfp;++m){
+ for(int n=0;n<p_Np;++n){
+ for(int f=0;f<p_Nfaces;++f){
+ f_LIFT[skL++] = mesh->LIFT[0][p_Nfp*p_Nfaces*n+(f+p_Nfaces*m)];
+ }
+ }
+ }
+#endif
+ cudaMalloc ((void**) &c_LIFT, sz);
+ cudaMemcpy( c_LIFT, f_LIFT, sz, cudaMemcpyHostToDevice);
+
+ /* Bind the array to the texture */
+ cudaBindTexture(0, t_LIFT, c_LIFT, sz);
+
+ /* DrDsDt */
+ sz = BSIZE*BSIZE*4*sizeof(float);
+
+ float* h_DrDsDt = (float*) calloc(BSIZE*BSIZE, sizeof(float4));
+ int sk = 0;
+ /* note transposed arrays to avoid "bank conflicts" */
+ for(int n=0;n<p_Np;++n){
+ for(int m=0;m<p_Np;++m){
+ h_DrDsDt[4*(m+n*BSIZE)+0] = mesh->Dr[0][n+m*p_Np];
+ h_DrDsDt[4*(m+n*BSIZE)+1] = mesh->Ds[0][n+m*p_Np];
+ h_DrDsDt[4*(m+n*BSIZE)+2] = mesh->Dt[0][n+m*p_Np];
+ }
+ }
+
+ cudaMalloc ((void**) &c_DrDsDt, sz);
+ cudaMemcpy( c_DrDsDt, h_DrDsDt, sz, cudaMemcpyHostToDevice);
+
+ /* Bind the array to the texture */
+ cudaBindTexture(0, t_DrDsDt, c_DrDsDt, sz);
+
+ free(h_DrDsDt);
+
+ /* vgeo */
+ double drdx, dsdx, dtdx;
+ double drdy, dsdy, dtdy;
+ double drdz, dsdz, dtdz, J;
+ float *vgeo = (float*) calloc(12*mesh->K, sizeof(float));
+
+ for(int k=0;k<mesh->K;++k){
+ GeometricFactors3d(mesh, k,
+ &drdx, &dsdx, &dtdx,
+ &drdy, &dsdy, &dtdy,
+ &drdz, &dsdz, &dtdz, &J);
+
+ vgeo[k*12+0] = drdx; vgeo[k*12+1] = drdy; vgeo[k*12+2] = drdz;
+ vgeo[k*12+4] = dsdx; vgeo[k*12+5] = dsdy; vgeo[k*12+6] = dsdz;
+ vgeo[k*12+8] = dtdx; vgeo[k*12+9] = dtdy; vgeo[k*12+10] = dtdz;
+
+ }
+
+ sz = mesh->K*12*sizeof(float);
+ cudaMalloc ((void**) &c_vgeo, sz);
+ cudaMemcpy( c_vgeo, vgeo, sz, cudaMemcpyHostToDevice);
+ cudaBindTexture(0, t_vgeo, c_vgeo, sz);
+
+ /* surfinfo (vmapM, vmapP, Fscale, Bscale, nx, ny, nz, 0) */
+ sz = mesh->K*p_Nfp*p_Nfaces*7*sizeof(float);
+ float* h_surfinfo = (float*) malloc(sz);
+
+ /* local-local info */
+ sk = 0;
+ int skP = -1;
+ double *nxk = BuildVector(mesh->Nfaces);
+ double *nyk = BuildVector(mesh->Nfaces);
+ double *nzk = BuildVector(mesh->Nfaces);
+ double *sJk = BuildVector(mesh->Nfaces);
+
+ double dt = 1e6;
+
+ for(int k=0;k<mesh->K;++k){
+
+ GeometricFactors3d(mesh, k,
+ &drdx, &dsdx, &dtdx,
+ &drdy, &dsdy, &dtdy,
+ &drdz, &dsdz, &dtdz, &J);
+
+ Normals3d(mesh, k, nxk, nyk, nzk, sJk);
+
+ for(int f=0;f<mesh->Nfaces;++f){
+
+ dt = min(dt, J/sJk[f]);
+
+ for(int m=0;m<p_Nfp;++m){
+ int n = m + f*p_Nfp + p_Nfp*p_Nfaces*k;
+ int idM = mesh->vmapM[n];
+ int idP = mesh->vmapP[n];
+ int nM = idM%p_Np;
+ int nP = idP%p_Np;
+ int kM = (idM-nM)/p_Np;
+ int kP = (idP-nP)/p_Np;
+ idM = nM + Nfields*BSIZE*kM;
+ idP = nP + Nfields*BSIZE*kP;
+
+ /* stub resolve some other way */
+ if(mesh->vmapP[n]<0){
+ idP = mesh->vmapP[n]; /* -ve numbers */
+ }
+
+ sk = 7*p_Nfp*p_Nfaces*k+m+f*p_Nfp;
+ h_surfinfo[sk + 0*p_Nfp*p_Nfaces] = idM;
+ h_surfinfo[sk + 1*p_Nfp*p_Nfaces] = idP;
+ h_surfinfo[sk + 2*p_Nfp*p_Nfaces] = sJk[f]/(2.*J);
+ h_surfinfo[sk + 3*p_Nfp*p_Nfaces] = (idM==idP)?-1.:1.;
+ h_surfinfo[sk + 4*p_Nfp*p_Nfaces] = nxk[f];
+ h_surfinfo[sk + 5*p_Nfp*p_Nfaces] = nyk[f];
+ h_surfinfo[sk + 6*p_Nfp*p_Nfaces] = nzk[f];
+ }
+ }
+ }
+
+ cudaMalloc ((void**) &c_surfinfo, sz);
+ cudaMemcpy( c_surfinfo, h_surfinfo, sz, cudaMemcpyHostToDevice);
+
+ cudaBindTexture(0, t_surfinfo, c_surfinfo, sz);
+
+ free(h_surfinfo);
+
+ sz = mesh->parNtotalout*sizeof(int);
+ cudaMalloc((void**) &(mesh->c_parmapOUT), sz);
+ cudaMemcpy(mesh->c_parmapOUT, mesh->parmapOUT, sz, cudaMemcpyHostToDevice);
+
+ return dt;
+}
+
+
+
+__global__ void MaxwellsGPU_VOL_Kernel3D(float *g_rhsQ){
+
+ /* fastest */
+ __device__ __shared__ float s_Q[p_Nfields*BSIZE];
+ __device__ __shared__ float s_facs[12];
+
+ const int n = threadIdx.x;
+ const int k = blockIdx.x;
+
+ /* "coalesced" */
+ int m = n+k*p_Nfields*BSIZE;
+ int id = n;
+ s_Q[id] = tex1Dfetch(t_Q, m); m+=BSIZE; id+=BSIZE;
+ s_Q[id] = tex1Dfetch(t_Q, m); m+=BSIZE; id+=BSIZE;
+ s_Q[id] = tex1Dfetch(t_Q, m); m+=BSIZE; id+=BSIZE;
+ s_Q[id] = tex1Dfetch(t_Q, m); m+=BSIZE; id+=BSIZE;
+ s_Q[id] = tex1Dfetch(t_Q, m); m+=BSIZE; id+=BSIZE;
+ s_Q[id] = tex1Dfetch(t_Q, m);
+
+ if(p_Np<12 && n==0)
+ for(m=0;m<12;++m)
+ s_facs[m] = tex1Dfetch(t_vgeo, 12*k+m);
+ else if(n<12 && p_Np>=12)
+ s_facs[n] = tex1Dfetch(t_vgeo, 12*k+n);
+
+ __syncthreads();
+
+ float dHxdr=0,dHxds=0,dHxdt=0;
+ float dHydr=0,dHyds=0,dHydt=0;
+ float dHzdr=0,dHzds=0,dHzdt=0;
+ float dExdr=0,dExds=0,dExdt=0;
+ float dEydr=0,dEyds=0,dEydt=0;
+ float dEzdr=0,dEzds=0,dEzdt=0;
+ float Q;
+
+ for(m=0;p_Np-m;){
+ float4 D = tex1Dfetch(t_DrDsDt, n+m*BSIZE);
+
+ id = m;
+ Q = s_Q[id]; dHxdr += D.x*Q; dHxds += D.y*Q; dHxdt += D.z*Q; id += BSIZE;
+ Q = s_Q[id]; dHydr += D.x*Q; dHyds += D.y*Q; dHydt += D.z*Q; id += BSIZE;
+ Q = s_Q[id]; dHzdr += D.x*Q; dHzds += D.y*Q; dHzdt += D.z*Q; id += BSIZE;
+ Q = s_Q[id]; dExdr += D.x*Q; dExds += D.y*Q; dExdt += D.z*Q; id += BSIZE;
+ Q = s_Q[id]; dEydr += D.x*Q; dEyds += D.y*Q; dEydt += D.z*Q; id += BSIZE;
+ Q = s_Q[id]; dEzdr += D.x*Q; dEzds += D.y*Q; dEzdt += D.z*Q;
+
+ ++m;
+#if ( (p_Np) % 2 )==0
+ D = tex1Dfetch(t_DrDsDt, n+m*BSIZE);
+
+ id = m;
+ Q = s_Q[id]; dHxdr += D.x*Q; dHxds += D.y*Q; dHxdt += D.z*Q; id += BSIZE;
+ Q = s_Q[id]; dHydr += D.x*Q; dHyds += D.y*Q; dHydt += D.z*Q; id += BSIZE;
+ Q = s_Q[id]; dHzdr += D.x*Q; dHzds += D.y*Q; dHzdt += D.z*Q; id += BSIZE;
+ Q = s_Q[id]; dExdr += D.x*Q; dExds += D.y*Q; dExdt += D.z*Q; id += BSIZE;
+ Q = s_Q[id]; dEydr += D.x*Q; dEyds += D.y*Q; dEydt += D.z*Q; id += BSIZE;
+ Q = s_Q[id]; dEzdr += D.x*Q; dEzds += D.y*Q; dEzdt += D.z*Q;
+
+ ++m;
+
+#if ( (p_Np)%3 )==0
+ D = tex1Dfetch(t_DrDsDt, n+m*BSIZE);
+
+ id = m;
+ Q = s_Q[id]; dHxdr += D.x*Q; dHxds += D.y*Q; dHxdt += D.z*Q; id += BSIZE;
+ Q = s_Q[id]; dHydr += D.x*Q; dHyds += D.y*Q; dHydt += D.z*Q; id += BSIZE;
+ Q = s_Q[id]; dHzdr += D.x*Q; dHzds += D.y*Q; dHzdt += D.z*Q; id += BSIZE;
+ Q = s_Q[id]; dExdr += D.x*Q; dExds += D.y*Q; dExdt += D.z*Q; id += BSIZE;
+ Q = s_Q[id]; dEydr += D.x*Q; dEyds += D.y*Q; dEydt += D.z*Q; id += BSIZE;
+ Q = s_Q[id]; dEzdr += D.x*Q; dEzds += D.y*Q; dEzdt += D.z*Q;
+
+ ++m;
+#endif
+#endif
+ }
+
+ const float drdx= s_facs[0];
+ const float drdy= s_facs[1];
+ const float drdz= s_facs[2];
+ const float dsdx= s_facs[4];
+ const float dsdy= s_facs[5];
+ const float dsdz= s_facs[6];
+ const float dtdx= s_facs[8];
+ const float dtdy= s_facs[9];
+ const float dtdz= s_facs[10];
+
+ m = n+p_Nfields*BSIZE*k;
+
+ g_rhsQ[m] = -(drdy*dEzdr+dsdy*dEzds+dtdy*dEzdt - drdz*dEydr-dsdz*dEyds-dtdz*dEydt); m += BSIZE;
+ g_rhsQ[m] = -(drdz*dExdr+dsdz*dExds+dtdz*dExdt - drdx*dEzdr-dsdx*dEzds-dtdx*dEzdt); m += BSIZE;
+ g_rhsQ[m] = -(drdx*dEydr+dsdx*dEyds+dtdx*dEydt - drdy*dExdr-dsdy*dExds-dtdy*dExdt); m += BSIZE;
+ g_rhsQ[m] = (drdy*dHzdr+dsdy*dHzds+dtdy*dHzdt - drdz*dHydr-dsdz*dHyds-dtdz*dHydt); m += BSIZE;
+ g_rhsQ[m] = (drdz*dHxdr+dsdz*dHxds+dtdz*dHxdt - drdx*dHzdr-dsdx*dHzds-dtdx*dHzdt); m += BSIZE;
+ g_rhsQ[m] = (drdx*dHydr+dsdx*dHyds+dtdx*dHydt - drdy*dHxdr-dsdy*dHxds-dtdy*dHxdt);
+}
+
+__global__ void MaxwellsGPU_SURF_Kernel3D(float *g_Q, float *g_rhsQ){
+
+ __device__ __shared__ float s_fluxQ[p_Nfields*p_Nfp*p_Nfaces];
+
+ const int n = threadIdx.x;
+ const int k = blockIdx.x;
+ int m;
+
+ /* grab surface nodes and store flux in shared memory */
+ if(n< (p_Nfp*p_Nfaces) ){
+ /* coalesced reads (maybe) */
+ m = 7*(k*p_Nfp*p_Nfaces)+n;
+ const int idM = tex1Dfetch(t_surfinfo, m); m += p_Nfp*p_Nfaces;
+ int idP = tex1Dfetch(t_surfinfo, m); m += p_Nfp*p_Nfaces;
+ const float Fsc = tex1Dfetch(t_surfinfo, m); m += p_Nfp*p_Nfaces;
+ const float Bsc = tex1Dfetch(t_surfinfo, m); m += p_Nfp*p_Nfaces;
+ const float nx = tex1Dfetch(t_surfinfo, m); m += p_Nfp*p_Nfaces;
+ const float ny = tex1Dfetch(t_surfinfo, m); m += p_Nfp*p_Nfaces;
+ const float nz = tex1Dfetch(t_surfinfo, m);
+
+ /* check if idP<0 */
+ double dHx, dHy, dHz, dEx, dEy, dEz;
+ if(idP<0){
+ idP = p_Nfields*(-1-idP);
+
+ dHx = Fsc*(tex1Dfetch(t_partQ, idP+0) - tex1Dfetch(t_Q, idM+0*BSIZE));
+ dHy = Fsc*(tex1Dfetch(t_partQ, idP+1) - tex1Dfetch(t_Q, idM+1*BSIZE));
+ dHz = Fsc*(tex1Dfetch(t_partQ, idP+2) - tex1Dfetch(t_Q, idM+2*BSIZE));
+
+ dEx = Fsc*(tex1Dfetch(t_partQ, idP+3) - tex1Dfetch(t_Q, idM+3*BSIZE));
+ dEy = Fsc*(tex1Dfetch(t_partQ, idP+4) - tex1Dfetch(t_Q, idM+4*BSIZE));
+ dEz = Fsc*(tex1Dfetch(t_partQ, idP+5) - tex1Dfetch(t_Q, idM+5*BSIZE));
+ }
+ else{
+ dHx = Fsc*(tex1Dfetch(t_Q, idP+0*BSIZE) - tex1Dfetch(t_Q, idM+0*BSIZE));
+ dHy = Fsc*(tex1Dfetch(t_Q, idP+1*BSIZE) - tex1Dfetch(t_Q, idM+1*BSIZE));
+ dHz = Fsc*(tex1Dfetch(t_Q, idP+2*BSIZE) - tex1Dfetch(t_Q, idM+2*BSIZE));
+
+ dEx = Fsc*(Bsc*tex1Dfetch(t_Q, idP+3*BSIZE) - tex1Dfetch(t_Q, idM+3*BSIZE));
+ dEy = Fsc*(Bsc*tex1Dfetch(t_Q, idP+4*BSIZE) - tex1Dfetch(t_Q, idM+4*BSIZE));
+ dEz = Fsc*(Bsc*tex1Dfetch(t_Q, idP+5*BSIZE) - tex1Dfetch(t_Q, idM+5*BSIZE));
+ }
+
+ const double ndotdH = nx*dHx + ny*dHy + nz*dHz;
+ const double ndotdE = nx*dEx + ny*dEy + nz*dEz;
+
+ m = n;
+ s_fluxQ[m] = -ny*dEz + nz*dEy + dHx - ndotdH*nx; m += p_Nfp*p_Nfaces;
+ s_fluxQ[m] = -nz*dEx + nx*dEz + dHy - ndotdH*ny; m += p_Nfp*p_Nfaces;
+ s_fluxQ[m] = -nx*dEy + ny*dEx + dHz - ndotdH*nz; m += p_Nfp*p_Nfaces;
+
+ s_fluxQ[m] = ny*dHz - nz*dHy + dEx - ndotdE*nx; m += p_Nfp*p_Nfaces;
+ s_fluxQ[m] = nz*dHx - nx*dHz + dEy - ndotdE*ny; m += p_Nfp*p_Nfaces;
+ s_fluxQ[m] = nx*dHy - ny*dHx + dEz - ndotdE*nz;
+ }
+
+ /* make sure all element data points are cached */
+ __syncthreads();
+
+ if(n< (p_Np))
+ {
+ float rhsHx = 0, rhsHy = 0, rhsHz = 0;
+ float rhsEx = 0, rhsEy = 0, rhsEz = 0;
+
+ int sk = n;
+ /* can manually unroll to 4 because there are 4 faces */
+ for(m=0;p_Nfaces*p_Nfp-m;){
+ const float4 L = tex1Dfetch(t_LIFT, sk); sk+=p_Np;
+
+ /* broadcast */
+ int sk1 = m;
+ rhsHx += L.x*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsHy += L.x*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsHz += L.x*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsEx += L.x*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsEy += L.x*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsEz += L.x*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ ++m;
+
+ /* broadcast */
+ sk1 = m;
+ rhsHx += L.y*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsHy += L.y*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsHz += L.y*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsEx += L.y*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsEy += L.y*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsEz += L.y*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ ++m;
+
+ /* broadcast */
+ sk1 = m;
+ rhsHx += L.z*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsHy += L.z*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsHz += L.z*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsEx += L.z*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsEy += L.z*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsEz += L.z*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ ++m;
+
+ /* broadcast */
+ sk1 = m;
+ rhsHx += L.w*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsHy += L.w*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsHz += L.w*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsEx += L.w*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsEy += L.w*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ rhsEz += L.w*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+ ++m;
+
+ }
+
+ m = n+p_Nfields*k*BSIZE;
+ g_rhsQ[m] += rhsHx; m += BSIZE;
+ g_rhsQ[m] += rhsHy; m += BSIZE;
+ g_rhsQ[m] += rhsHz; m += BSIZE;
+ g_rhsQ[m] += rhsEx; m += BSIZE;
+ g_rhsQ[m] += rhsEy; m += BSIZE;
+ g_rhsQ[m] += rhsEz; m += BSIZE;
+
+ }
+}
+
+
+__global__ void MaxwellsGPU_RK_Kernel3D(int Ntotal, float *g_resQ, float *g_rhsQ, float *g_Q, float fa, float fb, float fdt){
+
+ int n = blockIdx.x * blockDim.x + threadIdx.x;
+
+ if(n<Ntotal){
+ float rhs = g_rhsQ[n];
+ float res = g_resQ[n];
+ res = fa*res + fdt*rhs;
+
+ g_resQ[n] = res;
+ g_Q[n] += fb*res;
+ }
+
+}
+
+
+/* assumes data resides on device */
+void MaxwellsKernel3d(Mesh *mesh, float frka, float frkb, float fdt){
+
+ /* grab data from device and initiate sends */
+ MaxwellsMPISend3d(mesh);
+
+ int ThreadsPerBlock, BlocksPerGrid;
+
+ BlocksPerGrid = mesh->K;
+ ThreadsPerBlock = p_Np;
+
+ /* evaluate volume derivatives */
+ MaxwellsGPU_VOL_Kernel3D <<< BlocksPerGrid, ThreadsPerBlock >>> (c_rhsQ);
+
+ /* finalize sends and recvs, and transfer to device */
+ MaxwellsMPIRecv3d(mesh, c_partQ);
+
+ BlocksPerGrid = mesh->K;
+
+ if( ( p_Nfp*p_Nfaces ) > (p_Np) )
+ ThreadsPerBlock = p_Nfp*p_Nfaces;
+ else
+ ThreadsPerBlock = p_Np;
+
+ /* evaluate surface contributions */
+ MaxwellsGPU_SURF_Kernel3D <<< BlocksPerGrid, ThreadsPerBlock >>> (c_Q, c_rhsQ);
+
+ int Ntotal = mesh->K*BSIZE*p_Nfields;
+
+ ThreadsPerBlock = 256;
+ BlocksPerGrid = (Ntotal+ThreadsPerBlock-1)/ThreadsPerBlock;
+
+ /* update RK Step */
+ MaxwellsGPU_RK_Kernel3D<<< BlocksPerGrid, ThreadsPerBlock >>>
+ (Ntotal, c_resQ, c_rhsQ, c_Q, frka, frkb, fdt);
+
+}
+
+
+
+
+void gpu_set_data3d(int K,
+ double *d_Hx, double *d_Hy, double *d_Hz,
+ double *d_Ex, double *d_Ey, double *d_Ez){
+
+
+ float *f_Q = (float*) calloc(K*p_Nfields*BSIZE,sizeof(float));
+
+ /* also load into usual data matrices */
+
+ for(int k=0;k<K;++k){
+ int gk = k;
+ for(int n=0;n<p_Np;++n)
+ f_Q[n +k*BSIZE*p_Nfields] = d_Hx[n+gk*p_Np];
+ for(int n=0;n<p_Np;++n)
+ f_Q[n +BSIZE+k*BSIZE*p_Nfields] = d_Hy[n+gk*p_Np];
+ for(int n=0;n<p_Np;++n)
+ f_Q[n+2*BSIZE+k*BSIZE*p_Nfields] = d_Hz[n+gk*p_Np];
+ for(int n=0;n<p_Np;++n)
+ f_Q[n+3*BSIZE+k*BSIZE*p_Nfields] = d_Ex[n+gk*p_Np];
+ for(int n=0;n<p_Np;++n)
+ f_Q[n+4*BSIZE+k*BSIZE*p_Nfields] = d_Ey[n+gk*p_Np];
+ for(int n=0;n<p_Np;++n)
+ f_Q[n+5*BSIZE+k*BSIZE*p_Nfields] = d_Ez[n+gk*p_Np];
+ }
+
+ cudaMemcpy(c_Q, f_Q, BSIZE*K*p_Nfields*sizeof(float), cudaMemcpyHostToDevice);
+
+ free(f_Q);
+}
+
+void gpu_get_data3d(int K,
+ double *d_Hx, double *d_Hy, double *d_Hz,
+ double *d_Ex, double *d_Ey, double *d_Ez){
+
+ float *f_Q = (float*) calloc(K*p_Nfields*BSIZE,sizeof(float));
+
+ cudaMemcpy(f_Q, c_Q, K*BSIZE*p_Nfields*sizeof(float), cudaMemcpyDeviceToHost);
+
+ /* also load into usual data matrices */
+
+ for(int k=0;k<K;++k){
+ int gk = k;
+ for(int n=0;n<p_Np;++n)
+ d_Hx[n+gk*p_Np] = f_Q[n +k*BSIZE*p_Nfields];
+ for(int n=0;n<p_Np;++n)
+ d_Hy[n+gk*p_Np] = f_Q[n +BSIZE+k*BSIZE*p_Nfields];
+ for(int n=0;n<p_Np;++n)
+ d_Hz[n+gk*p_Np] = f_Q[n+2*BSIZE+k*BSIZE*p_Nfields];
+ for(int n=0;n<p_Np;++n)
+ d_Ex[n+gk*p_Np] = f_Q[n+3*BSIZE+k*BSIZE*p_Nfields];
+ for(int n=0;n<p_Np;++n)
+ d_Ey[n+gk*p_Np] = f_Q[n+4*BSIZE+k*BSIZE*p_Nfields];
+ for(int n=0;n<p_Np;++n)
+ d_Ez[n+gk*p_Np] = f_Q[n+5*BSIZE+k*BSIZE*p_Nfields];
+
+ }
+
+ free(f_Q);
+}
+
+__global__ void partial_get_kernel3d(int Ntotal, int *g_index, float *g_partQ){
+
+ int n = blockIdx.x * blockDim.x + threadIdx.x;
+
+ if(n<Ntotal)
+ g_partQ[n] = tex1Dfetch(t_Q, g_index[n]);
+
+}
+
+void get_partial_gpu_data3d(int Ntotal, int *g_index, float *h_partQ){
+
+ int ThreadsPerBlock = 256;
+ int BlocksPerGrid = (Ntotal+ThreadsPerBlock-1)/ThreadsPerBlock;
+
+ partial_get_kernel3d <<< BlocksPerGrid, ThreadsPerBlock >>> (Ntotal, g_index, c_tmp);
+
+ cudaMemcpy(h_partQ, c_tmp, Ntotal*sizeof(float), cudaMemcpyDeviceToHost);
+}
+
+}
diff --git a/benchmarks/CUDA/DG/src/MaxwellsMPI3d.c b/benchmarks/CUDA/DG/src/MaxwellsMPI3d.c
new file mode 100644
index 0000000..6bdf46f
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/MaxwellsMPI3d.c
@@ -0,0 +1,71 @@
+
+#include "mpi.h"
+#include <parmetisbin.h>
+
+#include "fem.h"
+
+static MPI_Request *mpi_out_requests = NULL;
+static MPI_Request *mpi_in_requests = NULL;
+
+static int Nmess = 0;
+
+void MaxwellsMPISend3d(Mesh *mesh){
+
+ int p;
+
+ int procid = mesh->procid;
+ int nprocs = mesh->nprocs;
+
+ MPI_Status status;
+
+ if(mpi_out_requests==NULL){
+ mpi_out_requests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ mpi_in_requests = (MPI_Request*) calloc(nprocs, sizeof(MPI_Request));
+ }
+
+#ifdef CUDA
+ get_partial_gpu_data3d(mesh->parNtotalout, mesh->c_parmapOUT, mesh->f_outQ);
+#endif
+
+ /* non-blocked send/recv partition surface data */
+ Nmess = 0;
+
+ /* now send piece to each proc */
+ int sk = 0;
+ for(p=0;p<nprocs;++p){
+
+ if(p!=procid){
+ int Nout = mesh->Npar[p]*p_Nfields*p_Nfp;
+ if(Nout){
+ /* symmetric communications (different ordering) */
+ MPI_Isend(mesh->f_outQ+sk, Nout, MPI_FLOAT, p, 6666+p, MPI_COMM_WORLD, mpi_out_requests +Nmess);
+ MPI_Irecv(mesh->f_inQ+sk, Nout, MPI_FLOAT, p, 6666+procid, MPI_COMM_WORLD, mpi_in_requests +Nmess);
+ sk+=Nout;
+ ++Nmess;
+ }
+ }
+ }
+
+}
+
+
+void MaxwellsMPIRecv3d(Mesh *mesh, float *c_partQ){
+ int p, n;
+ int nprocs = mesh->nprocs;
+
+ MPI_Status *instatus = (MPI_Status*) calloc(nprocs, sizeof(MPI_Status));
+ MPI_Status *outstatus = (MPI_Status*) calloc(nprocs, sizeof(MPI_Status));
+
+ MPI_Waitall(Nmess, mpi_in_requests, instatus);
+
+#ifdef CUDA
+ cudaMemcpy(c_partQ, mesh->f_inQ, mesh->parNtotalout*sizeof(float), cudaMemcpyHostToDevice);
+#endif
+
+ MPI_Waitall(Nmess, mpi_out_requests, outstatus);
+
+ free(outstatus);
+ free(instatus);
+
+}
+
diff --git a/benchmarks/CUDA/DG/src/MaxwellsRHS3d.c b/benchmarks/CUDA/DG/src/MaxwellsRHS3d.c
new file mode 100644
index 0000000..13ed5bf
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/MaxwellsRHS3d.c
@@ -0,0 +1,220 @@
+#include "mpi.h"
+#include "fem.h"
+
+void MaxwellsRHS3d(Mesh *mesh, float frka, float frkb, float fdt){
+
+ /* registers and temporary */
+ register unsigned int k, n;
+
+ /* mesh parameters */
+ const int K = mesh->K;
+
+ float *vgeo = mesh->vgeo;
+ float *surfinfo = mesh->surfinfo;
+ float *f_Dr = mesh->f_Dr;
+ float *f_Ds = mesh->f_Ds;
+ float *f_Dt = mesh->f_Dt;
+ float *f_LIFT = mesh->f_LIFT;
+
+ float *f_Q = mesh->f_Q;
+ float *f_rhsQ = mesh->f_rhsQ;
+ float *f_resQ = mesh->f_resQ;
+
+ float *f_inQ = mesh->f_inQ;
+ float *f_outQ = mesh->f_outQ;
+
+ int p;
+
+ /* mpi request buffer */
+ MPI_Request *mpi_out_requests = (MPI_Request*) calloc(mesh->nprocs, sizeof(MPI_Request));
+ MPI_Request *mpi_in_requests = (MPI_Request*) calloc(mesh->nprocs, sizeof(MPI_Request));
+
+ /* buffer outgoing node data */
+ for(n=0;n<mesh->parNtotalout;++n)
+ mesh->f_outQ[n] = f_Q[mesh->parmapOUT[n]];
+
+ /* do sends */
+ int sk = 0, Nmess = 0;
+ for(p=0;p<mesh->nprocs;++p){
+ if(p!=mesh->procid){
+ int Nout = mesh->Npar[p]*p_Nfields*p_Nfp;
+ if(Nout){
+ /* symmetric communications (different ordering) */
+ MPI_Isend(f_outQ+sk, Nout, MPI_FLOAT, p, 6666+p, MPI_COMM_WORLD, mpi_out_requests +Nmess);
+ MPI_Irecv(f_inQ+sk, Nout, MPI_FLOAT, p, 6666+mesh->procid, MPI_COMM_WORLD, mpi_in_requests +Nmess);
+ sk+=Nout;
+ ++Nmess;
+ }
+ }
+ }
+
+ for(k=0;k<K;++k){
+
+ /* NOTE: once k is known, all other indexing variables etc are derived */
+ register unsigned int n, m;
+
+ /* NOTE: should be local memory */
+ float Qk[p_Np*p_Nfields];
+
+ /* NOTE: index into geometric factors */
+ int geoid=k*12;
+
+ const float drdx = vgeo[geoid++], drdy = vgeo[geoid++], drdz = vgeo[geoid++]; ++geoid;
+ const float dsdx = vgeo[geoid++], dsdy = vgeo[geoid++], dsdz = vgeo[geoid++]; ++geoid;
+ const float dtdx = vgeo[geoid++], dtdy = vgeo[geoid++], dtdz = vgeo[geoid++]; ++geoid;
+
+ int id = k*p_Nfp*p_Nfaces;
+
+ /* NOTE: buffer element k into local storage */
+ float *qpt = f_Q+p_Nfields*p_Np*k;
+ for(m=0;m<p_Nfields*p_Np;++m){
+ Qk[m] = qpt[m];
+ }
+
+ int sk = 0;
+
+ for(n=0;n<p_Np;++n){
+
+ const float *ptDr = f_Dr+n*p_Np;
+ const float *ptDs = f_Ds+n*p_Np;
+ const float *ptDt = f_Dt+n*p_Np;
+
+ int skV = 0;
+
+ float rhsHx = 0, rhsHy = 0, rhsHz = 0;
+ float rhsEx = 0, rhsEy = 0, rhsEz = 0;
+
+ int sk3 = 0;
+ for(m=0;m<p_Np;++m){
+ const float dr = ptDr[m], ds = ptDs[m], dt = ptDt[m];
+ const float dx = drdx*dr+dsdx*ds+dtdx*dt;
+ const float dy = drdy*dr+dsdy*ds+dtdy*dt;
+ const float dz = drdz*dr+dsdz*ds+dtdz*dt;
+ const float nHx = Qk[sk3++];
+ const float nHy = Qk[sk3++];
+ const float nHz = Qk[sk3++];
+ const float nEx = Qk[sk3++];
+ const float nEy = Qk[sk3++];
+ const float nEz = Qk[sk3++];
+
+ rhsHx += dz*nEy-dy*nEz;
+ rhsHy += dx*nEz-dz*nEx;
+ rhsHz += dy*nEx-dx*nEy;
+
+ rhsEx += dy*nHz-dz*nHy;
+ rhsEy += dz*nHx-dx*nHz;
+ rhsEz += dx*nHy-dy*nHx;
+ }
+
+ int skR = p_Nfields*(n+k*p_Np);
+ f_rhsQ[skR++] = rhsHx;
+ f_rhsQ[skR++] = rhsHy;
+ f_rhsQ[skR++] = rhsHz;
+ f_rhsQ[skR++] = rhsEx;
+ f_rhsQ[skR++] = rhsEy;
+ f_rhsQ[skR++] = rhsEz;
+ }
+ }
+
+ /* DO RECV */
+ MPI_Status *instatus = (MPI_Status*) calloc(mesh->nprocs, sizeof(MPI_Status));
+ MPI_Waitall(Nmess, mpi_in_requests, instatus);
+ free(instatus);
+
+ for(k=0;k<K;++k){
+
+ /* NOTE: once k is known, all other indexing variables etc are derived */
+ register unsigned int n, m;
+
+ /* NOTE: should be local memory */
+ float fluxQ[p_Nfaces*p_Nfp*p_Nfields];
+
+ /* NOTE: index into geometric factors */
+ int surfid=k*7*p_Nfp*p_Nfaces;
+
+ int id = k*p_Nfp*p_Nfaces;
+
+ int sk = 0;
+
+ for(m=0;m<p_Nfp*p_Nfaces;++m){
+ int idM = surfinfo[surfid++];
+ int idP = surfinfo[surfid++];
+ const float FSc = surfinfo[surfid++];
+ const float BSc = surfinfo[surfid++];
+ const float NXf = surfinfo[surfid++];
+ const float NYf = surfinfo[surfid++];
+ const float NZf = surfinfo[surfid++];
+
+ float dHx, dHy, dHz, dEx, dEy, dEz;
+ if(idP<0){
+ idP = p_Nfields*(-1-idP);
+ dHx = FSc*(f_inQ[idP++] -f_Q[idM++]);
+ dHy = FSc*(f_inQ[idP++] -f_Q[idM++]);
+ dHz = FSc*(f_inQ[idP++] -f_Q[idM++]);
+ dEx = FSc*(f_inQ[idP++] -f_Q[idM++]);
+ dEy = FSc*(f_inQ[idP++] -f_Q[idM++]);
+ dEz = FSc*(f_inQ[idP] -f_Q[idM]);
+ }else{
+ dHx = FSc*( f_Q[idP++] -f_Q[idM++]);
+ dHy = FSc*( f_Q[idP++] -f_Q[idM++]);
+ dHz = FSc*( f_Q[idP++] -f_Q[idM++]);
+ dEx = FSc*(BSc*f_Q[idP++] -f_Q[idM++]);
+ dEy = FSc*(BSc*f_Q[idP++] -f_Q[idM++]);
+ dEz = FSc*(BSc*f_Q[idP] -f_Q[idM]);
+ }
+
+ const float ndotdH = NXf*dHx + NYf*dHy + NZf*dHz;
+ const float ndotdE = NXf*dEx + NYf*dEy + NZf*dEz;
+
+ fluxQ[sk++] = -NYf*dEz + NZf*dEy + dHx - ndotdH*NXf;
+ fluxQ[sk++] = -NZf*dEx + NXf*dEz + dHy - ndotdH*NYf;
+ fluxQ[sk++] = -NXf*dEy + NYf*dEx + dHz - ndotdH*NZf;
+
+ fluxQ[sk++] = NYf*dHz - NZf*dHy + dEx - ndotdE*NXf;
+ fluxQ[sk++] = NZf*dHx - NXf*dHz + dEy - ndotdE*NYf;
+ fluxQ[sk++] = NXf*dHy - NYf*dHx + dEz - ndotdE*NZf;
+
+ }
+
+ for(n=0;n<p_Np;++n){
+
+ const float *ptLIFT = f_LIFT+n*p_Nfp*p_Nfaces;
+
+ float rhsHx = 0, rhsHy = 0, rhsHz = 0;
+ float rhsEx = 0, rhsEy = 0, rhsEz = 0;
+ int sk3 = 0;
+
+ /* NOTE: hint for compiler (nvcc) to unroll */
+ for(m=0;m<p_Nfp*p_Nfaces;++m){
+ const float L = ptLIFT[m];
+ rhsHx += L*fluxQ[sk3++];
+ rhsHy += L*fluxQ[sk3++];
+ rhsHz += L*fluxQ[sk3++];
+ rhsEx += L*fluxQ[sk3++];
+ rhsEy += L*fluxQ[sk3++];
+ rhsEz += L*fluxQ[sk3++];
+ }
+
+ int skR = p_Nfields*(n+k*p_Np);
+ f_rhsQ[skR++] += rhsHx;
+ f_rhsQ[skR++] += rhsHy;
+ f_rhsQ[skR++] += rhsHz;
+ f_rhsQ[skR++] += rhsEx;
+ f_rhsQ[skR++] += rhsEy;
+ f_rhsQ[skR++] += rhsEz;
+ }
+ }
+
+ for(n=0;n<K*p_Np*p_Nfields;++n){
+ f_resQ[n] = frka*f_resQ[n]+fdt*f_rhsQ[n];
+ f_Q[n] += frkb*f_resQ[n];
+ }
+
+ /* make sure all messages went out */
+ MPI_Status *outstatus = (MPI_Status*) calloc(mesh->nprocs, sizeof(MPI_Status));
+ MPI_Waitall(Nmess, mpi_out_requests, outstatus);
+ free(outstatus);
+
+
+}
+
diff --git a/benchmarks/CUDA/DG/src/MaxwellsRun3d.c b/benchmarks/CUDA/DG/src/MaxwellsRun3d.c
new file mode 100644
index 0000000..d4ee9f9
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/MaxwellsRun3d.c
@@ -0,0 +1,60 @@
+#include "mpi.h"
+#include "fem.h"
+
+void MaxwellsRun3d(Mesh *mesh, double FinalTime, double dt){
+
+ double time = 0;
+ int INTRK, tstep=0;
+
+ double mpitime0 = MPI_Wtime();
+ MPI_Barrier(MPI_COMM_WORLD);
+
+ dt = 0.001;
+
+ /* outer time step loop */
+ while (time<FinalTime){
+
+ /* adjust final step to end exactly at FinalTime */
+ if (time+dt > FinalTime) { dt = FinalTime-time; }
+
+ for (INTRK=1; INTRK<=5; ++INTRK) {
+
+ /* compute rhs of TM-mode MaxwellsGPU's equations */
+ const float fdt = dt;
+ const float fa = (float)mesh->rk4a[INTRK-1];
+ const float fb = (float)mesh->rk4b[INTRK-1];
+
+#ifdef CUDA
+ MaxwellsKernel3d(mesh, fa, fb, fdt);
+#else
+ MaxwellsRHS3d(mesh, fa, fb, fdt);
+#endif
+ }
+
+ time += dt; /* increment current time */
+ tstep++; /* increment timestep */
+ }
+
+#ifdef CUDA
+ cudaThreadSynchronize();
+#endif
+
+ double flopsV = p_Np*p_Np*36 + p_Np*66; /* V3 */
+ double flopsS = p_Nfp*p_Nfaces*(15 + 10 + 36) + p_Np*(p_Nfaces*p_Nfp*12 + 6);
+ double flopsR = p_Np*p_Nfields*4;
+
+ int Kloc = mesh->K;
+
+ double mpitime1 = MPI_Wtime();
+
+ double time_total = mpitime1-mpitime0;
+
+ MPI_Barrier(MPI_COMM_WORLD);
+ printf("%d %d %lg %lg proc, N, time taken, GFLOPS/s (GPU)\n",
+ mesh->procid,
+ p_N,
+ time_total,
+ 5*(tstep)*( (flopsV+flopsS+flopsR)*((double)Kloc/(1.e9*time_total))));
+
+
+}
diff --git a/benchmarks/CUDA/DG/src/Mesh3d.c b/benchmarks/CUDA/DG/src/Mesh3d.c
new file mode 100644
index 0000000..af76acb
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/Mesh3d.c
@@ -0,0 +1,195 @@
+#include "mpi.h"
+#include "fem.h"
+
+Mesh *ReadMesh3d(char *filename){
+
+ int n;
+
+ Mesh *mesh = (Mesh*) calloc(1, sizeof(Mesh));
+
+ char buf[BUFSIZ];
+
+ FILE *fp = fopen(filename, "r");
+
+ /* assume modified Gambit neutral format */
+ for(n=0;n<6;++n)
+ fgets(buf, BUFSIZ, fp);
+
+ fgets(buf, BUFSIZ, fp);
+ sscanf(buf, "%d %d \n", &(mesh->Nv), &(mesh->K));
+ mesh->Nverts = 4; /* assume tets */
+ mesh->Nedges = 6; /* assume tets */
+ mesh->Nfaces = 4; /* assume tets */
+
+ fgets(buf, BUFSIZ, fp);
+ fgets(buf, BUFSIZ, fp);
+
+ /* read vertex coordinates */
+ double *VX = BuildVector(mesh->Nv);
+ double *VY = BuildVector(mesh->Nv);
+ double *VZ = BuildVector(mesh->Nv);
+ for(n=0;n<mesh->Nv;++n){
+ fgets(buf, BUFSIZ, fp);
+ sscanf(buf, "%*d %lf %lf %lf", VX+n, VY+n, VZ+n);
+ }
+
+ /* decide on parition */
+ int procid, nprocs;
+
+ MPI_Comm_rank(MPI_COMM_WORLD, &procid);
+ MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+
+ mesh->procid = procid;
+ mesh->nprocs = nprocs;
+
+ /* assume this proc owns a block of elements */
+
+ int Klocal, Kstart;
+ int *Kprocs = (int*) calloc(nprocs, sizeof(int));
+ int p;
+
+ int **newEToV, *newKprocs;
+ double **newVX, **newVY;
+
+ Klocal = (int) ( (double)(mesh->K)/(double)nprocs );
+
+ for(p=0;p<nprocs-1;++p){
+ Kprocs[p] = Klocal;
+ }
+ Kprocs[p] = Klocal + mesh->K - nprocs*Klocal;
+
+
+ Kstart= 0;
+ for(p=0;p<procid;++p)
+ Kstart += Kprocs[p];
+
+ Klocal = Kprocs[procid];
+
+ /* read element to vertex connectivity */
+ fgets(buf, BUFSIZ, fp);
+ fgets(buf, BUFSIZ, fp);
+ mesh->EToV = BuildIntMatrix(Klocal, mesh->Nverts);
+ mesh->GX = BuildMatrix(Klocal, mesh->Nverts);
+ mesh->GY = BuildMatrix(Klocal, mesh->Nverts);
+ mesh->GZ = BuildMatrix(Klocal, mesh->Nverts);
+
+ int sk = 0, v;
+ for(n=0;n<mesh->K;++n){
+ fgets(buf, BUFSIZ, fp);
+ if(n>=Kstart && n<Kstart+Klocal){
+ sscanf(buf, "%*d %*d %*d %d %d %d %d",
+ mesh->EToV[sk]+0, mesh->EToV[sk]+1,
+ mesh->EToV[sk]+2, mesh->EToV[sk]+3);
+
+ /* correct to 0-index */
+ --(mesh->EToV[sk][0]);
+ --(mesh->EToV[sk][1]);
+ --(mesh->EToV[sk][2]);
+ --(mesh->EToV[sk][3]);
+
+ for(v=0;v<mesh->Nverts;++v){
+ mesh->GX[sk][v] = VX[mesh->EToV[sk][v]];
+ mesh->GY[sk][v] = VY[mesh->EToV[sk][v]];
+ mesh->GZ[sk][v] = VZ[mesh->EToV[sk][v]];
+ }
+
+ ++sk;
+ }
+ }
+ fgets(buf, BUFSIZ, fp);
+ fgets(buf, BUFSIZ, fp);
+
+ mesh->K = Klocal;
+
+ fclose(fp);
+
+ return mesh;
+
+}
+
+void PrintMesh ( Mesh *mesh ){
+ int n;
+ printf("Mesh data: \n");
+ printf("\n K = %d\n", mesh->K);
+ printf("\n Nv = %d\n", mesh->Nv);
+ printf("\n Nverts = %d\n", mesh->Nverts);
+ printf("\n Node coordinates = \n");
+ printf("\n Element to vertex connectivity = \n");
+ for(n=0;n<mesh->K;++n){
+ printf("%d: %d %d %d %d\n", n,
+ mesh->EToV[n][0], mesh->EToV[n][1],
+ mesh->EToV[n][2], mesh->EToV[n][3]);
+ }
+
+}
+
+void GeometricFactors3d(Mesh *mesh, int k,
+ double *drdx, double *dsdx, double *dtdx,
+ double *drdy, double *dsdy, double *dtdy,
+ double *drdz, double *dsdz, double *dtdz,
+ double *J){
+
+ double x1 = mesh->GX[k][0], y1 = mesh->GY[k][0], z1 = mesh->GZ[k][0];
+ double x2 = mesh->GX[k][1], y2 = mesh->GY[k][1], z2 = mesh->GZ[k][1];
+ double x3 = mesh->GX[k][2], y3 = mesh->GY[k][2], z3 = mesh->GZ[k][2];
+ double x4 = mesh->GX[k][3], y4 = mesh->GY[k][3], z4 = mesh->GZ[k][3];
+
+ /* compute geometric factors of the following afine map */
+ /* x = 0.5*( (-1-r-s-t)*x1 + (1+r)*x2 + (1+s)*x3 + (1+t)*x4) */
+ /* y = 0.5*( (-1-r-s-t)*y1 + (1+r)*y2 + (1+s)*y3 + (1+t)*y4) */
+ /* z = 0.5*( (-1-r-s-t)*z1 + (1+r)*z2 + (1+s)*z3 + (1+t)*z4) */
+
+ double dxdr = (x2-x1)/2, dxds = (x3-x1)/2, dxdt = (x4-x1)/2;
+ double dydr = (y2-y1)/2, dyds = (y3-y1)/2, dydt = (y4-y1)/2;
+ double dzdr = (z2-z1)/2, dzds = (z3-z1)/2, dzdt = (z4-z1)/2;
+
+ *J =
+ dxdr*(dyds*dzdt-dzds*dydt)
+ -dydr*(dxds*dzdt-dzds*dxdt)
+ +dzdr*(dxds*dydt-dyds*dxdt);
+
+ *drdx = (dyds*dzdt - dzds*dydt)/(*J);
+ *drdy = -(dxds*dzdt - dzds*dxdt)/(*J);
+ *drdz = (dxds*dydt - dyds*dxdt)/(*J);
+
+ *dsdx = -(dydr*dzdt - dzdr*dydt)/(*J);
+ *dsdy = (dxdr*dzdt - dzdr*dxdt)/(*J);
+ *dsdz = -(dxdr*dydt - dydr*dxdt)/(*J);
+
+ *dtdx = (dydr*dzds - dzdr*dyds)/(*J);
+ *dtdy = -(dxdr*dzds - dzdr*dxds)/(*J);
+ *dtdz = (dxdr*dyds - dydr*dxds)/(*J);
+
+ if(*J<1e-10)
+ printf("warning: J = %lg\n", *J);
+
+}
+
+void Normals3d(Mesh *mesh, int k,
+ double *nx, double *ny, double *nz, double *sJ){
+
+ int f;
+
+ double drdx, dsdx, dtdx;
+ double drdy, dsdy, dtdy;
+ double drdz, dsdz, dtdz;
+ double J;
+
+ GeometricFactors3d(mesh, k,
+ &drdx, &dsdx, &dtdx,
+ &drdy, &dsdy, &dtdy,
+ &drdz, &dsdz, &dtdz,
+ &J);
+
+ nx[0] = -dtdx; nx[1] = -dsdx; nx[2] = drdx + dsdx + dtdx; nx[3] = -drdx;
+ ny[0] = -dtdy; ny[1] = -dsdy; ny[2] = drdy + dsdy + dtdy; ny[3] = -drdy;
+ nz[0] = -dtdz; nz[1] = -dsdz; nz[2] = drdz + dsdz + dtdz; nz[3] = -drdz;
+
+ for(f=0;f<4;++f){
+ sJ[f] = sqrt(nx[f]*nx[f]+ny[f]*ny[f]+nz[f]*nz[f]);
+ nx[f] /= sJ[f];
+ ny[f] /= sJ[f];
+ nz[f] /= sJ[f];
+ sJ[f] *= J;
+ }
+}
diff --git a/benchmarks/CUDA/DG/src/ParallelPairs.c b/benchmarks/CUDA/DG/src/ParallelPairs.c
new file mode 100644
index 0000000..7474c0b
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/ParallelPairs.c
@@ -0,0 +1,124 @@
+#include <strings.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "mpi.h"
+
+#define max(a,b) ( (a>b)?a:b )
+#define min(a,b) ( (a<b)?a:b )
+
+void ParallelPairs(void *objs, int Nmyobjs, int sizeobj,
+ int (*numget)(const void *),
+ void (*numset)(const void *, int ),
+ int (*procget)(const void *),
+ void (*marry)(const void *, const void *),
+ int (*compare_objs)(const void *, const void *)){
+
+ char *myobjs = (char*) objs;
+
+ int n, p, sk, cnt, num;
+ int maxind = 0;
+
+ int procid, nprocs;
+ MPI_Comm_rank(MPI_COMM_WORLD, &procid);
+ MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+
+ /* local sort */
+ qsort(myobjs, Nmyobjs, sizeobj, compare_objs);
+
+ /* TW: homework replace from here --------> */
+ maxind = 0;
+ for(n=0;n<Nmyobjs;++n){
+ num = numget(myobjs+sizeobj*n);
+ maxind = max(maxind, num);
+ }
+
+ int globalmaxind;
+ MPI_Allreduce(&maxind, &globalmaxind, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
+
+ int binsize = ceil( (double)(globalmaxind)/(double)nprocs ) + 10;
+
+ int *outN = (int*) calloc(nprocs, sizeof(int));
+ int *inN = (int*) calloc(nprocs, sizeof(int));
+
+ int *cumoutN = (int*) calloc(nprocs, sizeof(int));
+ int *cuminN = (int*) calloc(nprocs, sizeof(int));
+
+ sk = 0;
+ /* count the number of objs in each bin */
+ int binup = binsize;
+ for(p=0;p<nprocs;++p){
+ while( numget(myobjs+(sk*sizeobj) ) <= binup ){
+ ++(outN[p]);
+ ++sk;
+ if(sk==Nmyobjs){
+ break;
+ }
+ }
+ binup += binsize;
+ if(sk==Nmyobjs){
+ break;
+ }
+ }
+ /* TW: <---------- replace to here */
+
+ /* communicate numbers to be sent to each bin */
+ MPI_Alltoall(outN, 1, MPI_INT,
+ inN, 1, MPI_INT,
+ MPI_COMM_WORLD);
+
+ /* build incoming buffer */
+ int Notherobjs = 0;
+ for(p=0;p<nprocs;++p)
+ Notherobjs += inN[p];
+
+ for(p=0;p<nprocs;++p){
+ outN[p] *= sizeobj/sizeof(int);
+ inN[p] *= sizeobj/sizeof(int);
+ }
+
+ for(p=1;p<nprocs;++p){
+ cumoutN[p] = cumoutN[p-1]+outN[p-1];
+ cuminN[p] = cuminN[p-1] + inN[p-1];
+ }
+
+ /* fill up bins of objects from cloud */
+ char *otherobjs = (char*) calloc(Notherobjs*sizeobj, sizeof(char));
+
+ MPI_Alltoallv(myobjs, outN, cumoutN, MPI_INT,
+ otherobjs, inN, cuminN, MPI_INT,
+ MPI_COMM_WORLD);
+
+ /* sort the bin */
+ qsort(otherobjs, Notherobjs, sizeobj, compare_objs);
+
+ /* number unique objs consecutively in each bin */
+ for(n=1;n<Notherobjs;++n){
+ /* match ? */
+ if(!compare_objs(otherobjs+ n*sizeobj,
+ otherobjs+(n-1)*sizeobj)){
+
+ marry(otherobjs+n*sizeobj, otherobjs+(n-1)*sizeobj);
+ }
+ }
+
+ char *outobjs = (char*) calloc(Notherobjs*sizeobj, sizeof(char));
+ sk = 0;
+ for(p=0;p<nprocs;++p)
+ for(n=0;n<Notherobjs;++n)
+ if(procget(otherobjs+n*sizeobj)==p){
+ memcpy(outobjs+sk*sizeobj, otherobjs+n*sizeobj, sizeobj);
+ ++sk;
+ }
+
+ /* send results out */
+ MPI_Alltoallv(outobjs, inN, cuminN, MPI_INT,
+ myobjs, outN, cumoutN, MPI_INT,
+ MPI_COMM_WORLD);
+
+
+ free(otherobjs);
+ free(outN); free(inN); free(cumoutN); free(cuminN);
+}
+
diff --git a/benchmarks/CUDA/DG/src/StartUp3d.c b/benchmarks/CUDA/DG/src/StartUp3d.c
new file mode 100644
index 0000000..99281e5
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/StartUp3d.c
@@ -0,0 +1,123 @@
+
+#include "fem.h"
+
+void StartUp3d(Mesh *mesh){
+
+#if p_N==1
+#include "data3dN01.h"
+#elif p_N==2
+#include "data3dN02.h"
+#elif p_N==3
+#include "data3dN03.h"
+#elif p_N==4
+#include "data3dN04.h"
+#elif p_N==5
+#include "data3dN05.h"
+#elif p_N==6
+#include "data3dN06.h"
+#elif p_N==7
+#include "data3dN07.h"
+#elif p_N==8
+#include "data3dN08.h"
+#elif p_N==9
+#include "data3dN09.h"
+#endif
+
+ int n, m, k;
+
+ /* load r, s, t */
+ mesh->r = BuildVector(p_Np);
+ mesh->s = BuildVector(p_Np);
+ mesh->t = BuildVector(p_Np);
+ for(n=0;n<p_Np;++n){
+ mesh->r[n] = p_r[n];
+ mesh->s[n] = p_s[n];
+ mesh->t[n] = p_t[n];
+ }
+
+ /* load Dr, Ds, Dt */
+ mesh->Dr = BuildMatrix(p_Np, p_Np);
+ mesh->Ds = BuildMatrix(p_Np, p_Np);
+ mesh->Dt = BuildMatrix(p_Np, p_Np);
+ for(n=0;n<p_Np;++n){
+ for(m=0;m<p_Np;++m){
+ mesh->Dr[n][m] = p_Dr[n][m];
+ mesh->Ds[n][m] = p_Ds[n][m];
+ mesh->Dt[n][m] = p_Dt[n][m];
+ }
+ }
+
+ /* load LIFT */
+ mesh->LIFT = BuildMatrix(p_Np, p_Nfp*p_Nfaces);
+ for(n=0;n<p_Np;++n){
+ for(m=0;m<p_Nfp*p_Nfaces;++m){
+ mesh->LIFT[n][m] = p_LIFT[n][m];
+ }
+ }
+
+ mesh->Fmask = BuildIntMatrix(p_Nfaces, p_Nfp);
+ for(n=0;n<p_Nfaces;++n){
+ for(m=0;m<p_Nfp;++m){
+ mesh->Fmask[n][m] = p_Fmask[n][m];
+ }
+ }
+
+ /* low storage RK coefficients */
+ mesh->rk4a = BuildVector(5);
+ mesh->rk4a[0] = 0.0;
+ mesh->rk4a[1] = -567301805773.0 / 1357537059087.0;
+ mesh->rk4a[2] = -2404267990393.0 / 2016746695238.0;
+ mesh->rk4a[3] = -3550918686646.0 / 2091501179385.0;
+ mesh->rk4a[4] = -1275806237668.0 / 842570457699.0;
+
+ mesh->rk4b = BuildVector(5);
+ mesh->rk4b[0] = 1432997174477.0 / 9575080441755.0;
+ mesh->rk4b[1] = 5161836677717.0 / 13612068292357.0;
+ mesh->rk4b[2] = 1720146321549.0 / 2090206949498.0;
+ mesh->rk4b[3] = 3134564353537.0 / 4481467310338.0;
+ mesh->rk4b[4] = 2277821191437.0 / 14882151754819.0;
+
+ mesh->rk4c = BuildVector(6);
+ mesh->rk4c[0] = 0.0;
+ mesh->rk4c[1] = 1432997174477.0 / 9575080441755.0;
+ mesh->rk4c[2] = 2526269341429.0 / 6820363962896.0;
+ mesh->rk4c[3] = 2006345519317.0 / 3224310063776.0;
+ mesh->rk4c[4] = 2802321613138.0 / 2924317926251.0;
+ mesh->rk4c[5] = 1.0;
+
+ /* build coordinates */
+ mesh->x = BuildMatrix(mesh->K, p_Np);
+ mesh->y = BuildMatrix(mesh->K, p_Np);
+ mesh->z = BuildMatrix(mesh->K, p_Np);
+
+ for(k=0;k<mesh->K;++k){
+ for(n=0;n<p_Np;++n){
+ double r = mesh->r[n];
+ double s = mesh->s[n];
+ double t = mesh->t[n];
+
+ mesh->x[k][n] = 0.5*( -mesh->GX[k][0]*(r+s+t+1) +
+ mesh->GX[k][1]*(1.+r) +
+ mesh->GX[k][2]*(1.+s) +
+ mesh->GX[k][3]*(1.+t)
+ );
+
+ mesh->y[k][n] = 0.5*( -mesh->GY[k][0]*(r+s+t+1) +
+ mesh->GY[k][1]*(1.+r) +
+ mesh->GY[k][2]*(1.+s) +
+ mesh->GY[k][3]*(1.+t)
+ );
+
+ mesh->z[k][n] = 0.5*( -mesh->GZ[k][0]*(r+s+t+1) +
+ mesh->GZ[k][1]*(1.+r) +
+ mesh->GZ[k][2]*(1.+s) +
+ mesh->GZ[k][3]*(1.+t)
+ );
+ }
+ }
+
+ /* build node-node connectivity maps */
+ void BuildMaps3d(Mesh *mesh);
+ BuildMaps3d(mesh);
+
+}
diff --git a/benchmarks/CUDA/DG/src/Utils.c b/benchmarks/CUDA/DG/src/Utils.c
new file mode 100644
index 0000000..e0364cf
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/Utils.c
@@ -0,0 +1,155 @@
+#include "fem.h"
+
+/* some very basic memory allocation routines */
+
+/* row major storage for a 2D matrix array */
+double **BuildMatrix(int Nrows, int Ncols){
+ int n;
+ double **A = (double**) calloc(Nrows, sizeof(double*));
+
+ A[0] = (double*) calloc(Nrows*Ncols, sizeof(double));
+
+ for(n=1;n<Nrows;++n){
+ A[n] = A[n-1]+ Ncols;
+ }
+
+ return A;
+}
+
+double *BuildVector(int Nrows){
+
+ double *A = (double*) calloc(Nrows, sizeof(double));
+
+ return A;
+}
+
+/* row major storage for a 2D matrix array */
+int **BuildIntMatrix(int Nrows, int Ncols){
+ int n;
+ int **A = (int**) calloc(Nrows, sizeof(int*));
+
+ A[0] = (int*) calloc(Nrows*Ncols, sizeof(int));
+
+ for(n=1;n<Nrows;++n){
+ A[n] = A[n-1]+ Ncols;
+ }
+
+ return A;
+}
+
+int *BuildIntVector(int Nrows){
+
+ int *A = (int*) calloc(Nrows, sizeof(int));
+
+ return A;
+}
+
+double *DestroyVector(double *v){
+ free(v);
+ return NULL;
+}
+
+double **DestroyMatrix(double **A){
+ free(A[0]);
+ free(A);
+
+ return NULL;
+}
+
+int *DestroyIntVector(int *v){
+ free(v);
+ return NULL;
+}
+
+int **DestroyIntMatrix(int **A){
+ free(A[0]);
+ free(A);
+
+ return NULL;
+}
+
+void PrintMatrix(char *message, double **A, int Nrows, int Ncols){
+ int n,m;
+
+ printf("%s\n", message);
+ for(n=0;n<Nrows;++n){
+ for(m=0;m<Ncols;++m){
+ printf(" %g ", A[n][m]);
+ }
+ printf(" \n");
+ }
+}
+
+
+void SaveMatrix(char *filename, double **A, int Nrows, int Ncols){
+ int n,m;
+
+ FILE *fp = fopen(filename, "w");
+
+ for(n=0;n<Nrows;++n){
+ for(m=0;m<Ncols;++m){
+ fprintf(fp, " %g ", A[n][m]);
+ }
+ fprintf(fp, " \n");
+ }
+
+ fclose(fp);
+}
+
+
+int trianglebase(Mesh *mesh, int k){
+
+ double x1 = mesh->GX[k][0];
+ double x2 = mesh->GX[k][1];
+ double x3 = mesh->GX[k][2];
+
+ double y1 = mesh->GY[k][0];
+ double y2 = mesh->GY[k][1];
+ double y3 = mesh->GY[k][2];
+
+ double d1 = (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2);
+ double d2 = (x2-x3)*(x2-x3) + (y2-y3)*(y2-y3);
+ double d3 = (x3-x1)*(x3-x1) + (y3-y1)*(y3-y1);
+
+ /* find maximum length face */
+ if(d1>=d2 && d1>=d3)
+ return 0;
+ else if(d2>=d3)
+ return 1;
+
+ return 2;
+}
+
+
+int tetbase(Mesh *mesh, int k){
+
+ double x1 = mesh->GX[k][0];
+ double x2 = mesh->GX[k][1];
+ double x3 = mesh->GX[k][2];
+ double x4 = mesh->GX[k][3];
+
+ double y1 = mesh->GY[k][0];
+ double y2 = mesh->GY[k][1];
+ double y3 = mesh->GY[k][2];
+ double y4 = mesh->GY[k][3];
+
+ double z1 = mesh->GZ[k][0];
+ double z2 = mesh->GZ[k][1];
+ double z3 = mesh->GZ[k][2];
+ double z4 = mesh->GZ[k][3];
+
+ double d1 = (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2);
+ double d2 = (x2-x3)*(x2-x3) + (y2-y3)*(y2-y3) + (z2-z3)*(z2-z3);
+ double d3 = (x3-x4)*(x3-x4) + (y3-y4)*(y3-y4) + (z3-z4)*(z3-z4);
+ double d4 = (x4-x1)*(x4-x1) + (y4-y1)*(y4-y1) + (z4-z1)*(z4-z1);
+
+ /* find maximum length face */
+ if(d1>=d2 && d1>=d3 && d1>=d4)
+ return 0;
+ else if(d2>=d3 && d2>=d4)
+ return 1;
+ else if(d3>=d4)
+ return 2;
+
+ return 3;
+}
diff --git a/benchmarks/CUDA/DG/src/buildsource.m b/benchmarks/CUDA/DG/src/buildsource.m
new file mode 100644
index 0000000..c93d363
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/buildsource.m
@@ -0,0 +1,106 @@
+
+function buildsource(in_N)
+
+Globals2D;
+
+
+NODETOL = 1e-8;
+
+N = in_N;
+
+Np = (N+1)*(N+2)/2; Nfp = N+1; Nfaces = 3;
+
+[x,y] = Nodes2D(N);
+[r,s] = xytors(x,y);
+
+% find all the nodes that lie on each edge
+fmask1 = find( abs(s+1) < NODETOL)';
+fmask2 = find( abs(r+s) < NODETOL)';
+fmask3 = find( abs(r+1) < NODETOL)';
+Fmask = [fmask1;fmask2;fmask3]';
+
+% vandermonde matrix
+V = Vandermonde2D(N, r, s);
+
+% derivative matrices
+[Dr, Ds] = Dmatrices2D(N, r, s, V);
+
+% lift matrix
+LIFT = Lift2D();
+
+% correct to 0-index
+Fmask = Fmask-1;
+
+
+fid = fopen(sprintf('dataN%02d.h', N), 'w');
+
+fprintf(fid, sprintf('#ifndef DATAN%02d \n\n', N));
+fprintf(fid, sprintf('#define DATAN%02d 1\n\n', N));
+fprintf(fid, 'double p_r[%d] = { ', Np);
+fprintf(fid, '%17.15g ', r(1));
+for n=2:Np
+ fprintf(fid, ', %17.15g ', r(n));
+end
+fprintf(fid, '};\n')
+
+fprintf(fid, 'double p_s[%d] = {', Np);
+fprintf(fid, '%17.15g ', s(1));
+for n=2:Np
+ fprintf(fid, ', %17.15g ', s(n));
+end
+fprintf(fid, '};\n')
+
+fprintf(fid, 'double p_Dr[%d][%d] = {', Np, Np);
+for m=1:Np
+fprintf(fid, '{%17.15g ', Dr(m,1));
+for n=2:Np
+ fprintf(fid, ', %17.15g ', Dr(m,n));
+end
+if(m<Np)
+fprintf(fid, '},\n')
+else
+fprintf(fid, '}};\n')
+end
+end
+
+fprintf(fid, 'double p_Ds[%d][%d] = {', Np, Np);
+for m=1:Np
+fprintf(fid, '{%17.15g ', Ds(m,1));
+for n=2:Np
+ fprintf(fid, ', %17.15g ', Ds(m,n));
+end
+if(m<Np)
+fprintf(fid, '},\n')
+else
+fprintf(fid, '}};\n')
+end
+end
+
+fprintf(fid, 'double p_LIFT[%d][%d] = {', Np, Nfp*Nfaces);
+for m=1:Np
+fprintf(fid, '{%17.15g ', LIFT(m,1));
+for n=2:Nfp*Nfaces
+ fprintf(fid, ', %17.15g ', LIFT(m,n));
+end
+if(m<Np)
+fprintf(fid, '},\n')
+else
+fprintf(fid, '}};\n')
+end
+end
+
+fprintf(fid, 'int p_Fmask[%d][%d] = {', Nfaces, Nfp);
+for m=1:Nfaces
+fprintf(fid, '{%d ', Fmask(1,m));
+for n=2:Nfp
+ fprintf(fid, ', %d ', Fmask(n,m));
+end
+if(m<Nfaces)
+fprintf(fid, '},\n')
+else
+fprintf(fid, '}};\n')
+end
+end
+
+fprintf(fid, sprintf('\n#endif\n'));
+fclose(fid)
diff --git a/benchmarks/CUDA/DG/src/gpuDG.vpj b/benchmarks/CUDA/DG/src/gpuDG.vpj
new file mode 100644
index 0000000..0142803
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/gpuDG.vpj
@@ -0,0 +1,173 @@
+<!DOCTYPE Project SYSTEM "http://www.slickedit.com/dtd/vse/10.0/vpj.dtd">
+<Project
+ Version="10.0"
+ VendorName="SlickEdit"
+ WorkingDir=".">
+ <Config
+ Name="Debug"
+ Type="cpp"
+ DebugCallbackName="gdb"
+ OutputFile="%bdgpuDG"
+ CompilerConfigName="Latest Version">
+ <Menu>
+ <Target
+ Name="Compile"
+ MenuCaption="&amp;Compile"
+ CaptureOutputWith="ProcessBuffer"
+ OutputExts="*.o"
+ SaveOption="SaveCurrent"
+ RunFromDir="%rw">
+ <Exec CmdLine='cc -c -g %i %defd -o "%bd%n.o" "%f"'/>
+ </Target>
+ <Target
+ Name="Link"
+ MenuCaption="&amp;Link">
+ <Exec CmdLine='cc -g -o "%o" %f %libs'/>
+ </Target>
+ <Target
+ Name="Build"
+ MenuCaption="&amp;Build"
+ CaptureOutputWith="ProcessBuffer"
+ SaveOption="SaveWorkspaceFiles"
+ RunFromDir="%rw">
+ <Exec CmdLine='"%(VSLICKBIN1)vsbuild" build "%w" "%r"'/>
+ </Target>
+ <Target
+ Name="Rebuild"
+ MenuCaption="&amp;Rebuild"
+ CaptureOutputWith="ProcessBuffer"
+ SaveOption="SaveWorkspaceFiles"
+ RunFromDir="%rw">
+ <Exec CmdLine='"%(VSLICKBIN1)vsbuild" rebuild "%w" "%r"'/>
+ </Target>
+ <Target
+ Name="Debug"
+ MenuCaption="&amp;Debug"
+ SaveOption="SaveNone"
+ BuildFirst="1"
+ CaptureOutputWith="ProcessBuffer"
+ RunFromDir="%rw">
+ <Exec CmdLine='vsdebugio -prog "%o"'/>
+ </Target>
+ <Target
+ Name="Execute"
+ MenuCaption="E&amp;xecute"
+ SaveOption="SaveNone"
+ BuildFirst="1"
+ CaptureOutputWith="ProcessBuffer"
+ RunFromDir="%rw">
+ <Exec CmdLine='"%o"'/>
+ </Target>
+ </Menu>
+ </Config>
+ <Config
+ Name="Release"
+ Type="cpp"
+ DebugCallbackName="gdb"
+ OutputFile="%bdgpuDG"
+ CompilerConfigName="Latest Version">
+ <Menu>
+ <Target
+ Name="Compile"
+ MenuCaption="&amp;Compile"
+ CaptureOutputWith="ProcessBuffer"
+ SaveOption="SaveCurrent"
+ OutputExts="*.o"
+ RunFromDir="%rw">
+ <Exec CmdLine='cc -c -O %i %defd -o "%bd%n.o" "%f"'/>
+ </Target>
+ <Target
+ Name="Link"
+ MenuCaption="&amp;Link">
+ <Exec CmdLine='cc -O -o "%o" %f %libs'/>
+ </Target>
+ <Target
+ Name="Build"
+ MenuCaption="&amp;Build"
+ CaptureOutputWith="ProcessBuffer"
+ SaveOption="SaveWorkspaceFiles"
+ RunFromDir="%rw">
+ <Exec CmdLine='"%(VSLICKBIN1)vsbuild" build "%w" "%r"'/>
+ </Target>
+ <Target
+ Name="Rebuild"
+ MenuCaption="&amp;Rebuild"
+ CaptureOutputWith="ProcessBuffer"
+ SaveOption="SaveWorkspaceFiles"
+ RunFromDir="%rw">
+ <Exec CmdLine='"%(VSLICKBIN1)vsbuild" rebuild "%w" "%r"'/>
+ </Target>
+ <Target
+ Name="Debug"
+ MenuCaption="&amp;Debug"
+ SaveOption="SaveNone"
+ BuildFirst="1"
+ CaptureOutputWith="ProcessBuffer"
+ RunFromDir="%rw">
+ <Exec CmdLine='vsdebugio -prog "%o"'/>
+ </Target>
+ <Target
+ Name="Execute"
+ MenuCaption="E&amp;xecute"
+ SaveOption="SaveNone"
+ BuildFirst="1"
+ CaptureOutputWith="ProcessBuffer"
+ RunFromDir="%rw">
+ <Exec CmdLine='"%o"'/>
+ </Target>
+ </Menu>
+ </Config>
+ <Rules Name="Compile">
+ <Rule
+ InputExts="*.s"
+ OutputExts="*.o"
+ LinkObject="1">
+ <Exec CmdLine='as -o "%bd%n.o" "%f"'/>
+ </Rule>
+ </Rules>
+ <Files>
+ <Folder
+ Name="Source Files"
+ Filters="*.c;*.C;*.cc;*.cpp;*.cp;*.cxx;*.prg;*.pas;*.dpr;*.asm;*.s;*.bas;*.java;*.cs;*.sc;*.e;*.cob;*.html;*.rc;*.tcl;*.py;*.pl">
+ <F N="BuildMaps2d.c"/>
+ <F N="BuildMaps3d.c"/>
+ <F N="FacePair2d.c"/>
+ <F N="FacePair3d.c"/>
+ <F N="InitCPU2d.c"/>
+ <F N="InitCPU3d.c"/>
+ <F N="LoadBalance2d.c"/>
+ <F N="LoadBalance3d.c"/>
+ <F N="MaxwellsDriver2d.c"/>
+ <F N="MaxwellsDriver3d.c"/>
+ <F N="MaxwellsMPI2d.c"/>
+ <F N="MaxwellsMPI3d.c"/>
+ <F N="MaxwellsRHS2d.c"/>
+ <F N="MaxwellsRHS3d.c"/>
+ <F N="MaxwellsRun2d.c"/>
+ <F N="MaxwellsRun3d.c"/>
+ <F N="Mesh2d.c"/>
+ <F N="Mesh3d.c"/>
+ <F N="ParallelPairs.c"/>
+ <F N="StartUp2d.c"/>
+ <F N="StartUp3d.c"/>
+ <F N="Utils.c"/>
+ </Folder>
+ <Folder
+ Name="Header Files"
+ Filters="*.h;*.H;*.hh;*.hpp;*.hxx;*.inc;*.sh;*.cpy;*.if"/>
+ <Folder
+ Name="Resource Files"
+ Filters="*.ico;*.cur;*.dlg"/>
+ <Folder
+ Name="Bitmaps"
+ Filters="*.bmp"/>
+ <Folder
+ Name="Other Files"
+ Filters="">
+ <F N="buildsource.m"/>
+ <F N="MaxwellsKernel2d.cu"/>
+ <F N="MaxwellsKernel3d.cu"/>
+ <F N="tictoc.cu"/>
+ </Folder>
+ </Files>
+</Project>
diff --git a/benchmarks/CUDA/DG/src/tictoc.cu b/benchmarks/CUDA/DG/src/tictoc.cu
new file mode 100644
index 0000000..f1335c3
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/tictoc.cu
@@ -0,0 +1,23 @@
+cudaEvent_t cstart, cstop;
+
+void cudatic(){
+
+ cudaEventCreate(&cstart);
+ cudaEventCreate(&cstop);
+
+ cudaEventRecord(cstart, 0);
+}
+
+float cudatoc(){
+
+ cudaEventRecord(cstop, 0);
+ cudaEventSynchronize(cstop);
+ float elapsedTime;
+ cudaEventElapsedTime(&elapsedTime, cstart, cstop);
+
+ /* return elapsed time in seconds */
+ return elapsedTime/1000.0;
+
+}
+
+
generated by cgit v1.3 (git 2.53.0) at 2026-07-10 23:31:13 +0000