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, 0 insertions, 3171 deletions
diff --git a/benchmarks/CUDA/DG/src/BuildMaps2d.c b/benchmarks/CUDA/DG/src/BuildMaps2d.c
deleted file mode 100644
index 7d0d57c..0000000
--- a/benchmarks/CUDA/DG/src/BuildMaps2d.c
+++ /dev/null
@@ -1,231 +0,0 @@
-#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
deleted file mode 100644
index 89f7a44..0000000
--- a/benchmarks/CUDA/DG/src/BuildMaps3d.c
+++ /dev/null
@@ -1,279 +0,0 @@
-#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
deleted file mode 100644
index 92e47f5..0000000
--- a/benchmarks/CUDA/DG/src/FacePair2d.c
+++ /dev/null
@@ -1,154 +0,0 @@
-#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
deleted file mode 100644
index 41d1fad..0000000
--- a/benchmarks/CUDA/DG/src/FacePair3d.c
+++ /dev/null
@@ -1,201 +0,0 @@
-#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
deleted file mode 100644
index 2fe907f..0000000
--- a/benchmarks/CUDA/DG/src/InitCPU3d.c
+++ /dev/null
@@ -1,152 +0,0 @@
-#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
deleted file mode 100644
index 586cd8c..0000000
--- a/benchmarks/CUDA/DG/src/LoadBalance3d.c
+++ /dev/null
@@ -1,243 +0,0 @@
-#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
deleted file mode 100644
index dd761c7..0000000
--- a/benchmarks/CUDA/DG/src/MaxwellsDriver3d.c
+++ /dev/null
@@ -1,110 +0,0 @@
-#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
deleted file mode 100644
index bfa4396..0000000
--- a/benchmarks/CUDA/DG/src/MaxwellsKernel3d.cu
+++ /dev/null
@@ -1,551 +0,0 @@
-/* -*- 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
deleted file mode 100644
index 6bdf46f..0000000
--- a/benchmarks/CUDA/DG/src/MaxwellsMPI3d.c
+++ /dev/null
@@ -1,71 +0,0 @@
-
-#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
deleted file mode 100644
index 13ed5bf..0000000
--- a/benchmarks/CUDA/DG/src/MaxwellsRHS3d.c
+++ /dev/null
@@ -1,220 +0,0 @@
-#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
deleted file mode 100644
index d4ee9f9..0000000
--- a/benchmarks/CUDA/DG/src/MaxwellsRun3d.c
+++ /dev/null
@@ -1,60 +0,0 @@
-#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
deleted file mode 100644
index af76acb..0000000
--- a/benchmarks/CUDA/DG/src/Mesh3d.c
+++ /dev/null
@@ -1,195 +0,0 @@
-#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
deleted file mode 100644
index 7474c0b..0000000
--- a/benchmarks/CUDA/DG/src/ParallelPairs.c
+++ /dev/null
@@ -1,124 +0,0 @@
-#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
deleted file mode 100644
index 99281e5..0000000
--- a/benchmarks/CUDA/DG/src/StartUp3d.c
+++ /dev/null
@@ -1,123 +0,0 @@
-
-#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
deleted file mode 100644
index e0364cf..0000000
--- a/benchmarks/CUDA/DG/src/Utils.c
+++ /dev/null
@@ -1,155 +0,0 @@
-#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
deleted file mode 100644
index c93d363..0000000
--- a/benchmarks/CUDA/DG/src/buildsource.m
+++ /dev/null
@@ -1,106 +0,0 @@
-
-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
deleted file mode 100644
index 0142803..0000000
--- a/benchmarks/CUDA/DG/src/gpuDG.vpj
+++ /dev/null
@@ -1,173 +0,0 @@
-<!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
deleted file mode 100644
index f1335c3..0000000
--- a/benchmarks/CUDA/DG/src/tictoc.cu
+++ /dev/null
@@ -1,23 +0,0 @@
-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-11 14:56:17 +0000