From 11b308e7363e937966b035b4891db32b4eece3bf Mon Sep 17 00:00:00 2001 From: Tor Aamodt Date: Fri, 1 Oct 2010 08:55:28 -0800 Subject: integrating recent changes from fermi-test into fermi (i'll use "fermi" for more disruptive changes to the pipeline model such as updating the MSHRs and getting rid of the warp tracker, ripping out DWF, etc...) [git-p4: depot-paths = "//depot/gpgpu_sim_research/fermi/distribution/": change = 7805] --- benchmarks/CUDA/DG/src/BuildMaps2d.c | 231 ------------ benchmarks/CUDA/DG/src/BuildMaps3d.c | 279 --------------- benchmarks/CUDA/DG/src/FacePair2d.c | 154 -------- benchmarks/CUDA/DG/src/FacePair3d.c | 201 ----------- benchmarks/CUDA/DG/src/InitCPU3d.c | 152 -------- benchmarks/CUDA/DG/src/LoadBalance3d.c | 243 ------------- benchmarks/CUDA/DG/src/MaxwellsDriver3d.c | 110 ------ benchmarks/CUDA/DG/src/MaxwellsKernel3d.cu | 551 ----------------------------- benchmarks/CUDA/DG/src/MaxwellsMPI3d.c | 71 ---- benchmarks/CUDA/DG/src/MaxwellsRHS3d.c | 220 ------------ benchmarks/CUDA/DG/src/MaxwellsRun3d.c | 60 ---- benchmarks/CUDA/DG/src/Mesh3d.c | 195 ---------- benchmarks/CUDA/DG/src/ParallelPairs.c | 124 ------- benchmarks/CUDA/DG/src/StartUp3d.c | 123 ------- benchmarks/CUDA/DG/src/Utils.c | 155 -------- benchmarks/CUDA/DG/src/buildsource.m | 106 ------ benchmarks/CUDA/DG/src/gpuDG.vpj | 173 --------- benchmarks/CUDA/DG/src/tictoc.cu | 23 -- 18 files changed, 3171 deletions(-) delete mode 100644 benchmarks/CUDA/DG/src/BuildMaps2d.c delete mode 100644 benchmarks/CUDA/DG/src/BuildMaps3d.c delete mode 100644 benchmarks/CUDA/DG/src/FacePair2d.c delete mode 100644 benchmarks/CUDA/DG/src/FacePair3d.c delete mode 100644 benchmarks/CUDA/DG/src/InitCPU3d.c delete mode 100644 benchmarks/CUDA/DG/src/LoadBalance3d.c delete mode 100644 benchmarks/CUDA/DG/src/MaxwellsDriver3d.c delete mode 100644 benchmarks/CUDA/DG/src/MaxwellsKernel3d.cu delete mode 100644 benchmarks/CUDA/DG/src/MaxwellsMPI3d.c delete mode 100644 benchmarks/CUDA/DG/src/MaxwellsRHS3d.c delete mode 100644 benchmarks/CUDA/DG/src/MaxwellsRun3d.c delete mode 100644 benchmarks/CUDA/DG/src/Mesh3d.c delete mode 100644 benchmarks/CUDA/DG/src/ParallelPairs.c delete mode 100644 benchmarks/CUDA/DG/src/StartUp3d.c delete mode 100644 benchmarks/CUDA/DG/src/Utils.c delete mode 100644 benchmarks/CUDA/DG/src/buildsource.m delete mode 100644 benchmarks/CUDA/DG/src/gpuDG.vpj delete mode 100644 benchmarks/CUDA/DG/src/tictoc.cu (limited to 'benchmarks/CUDA/DG/src') 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 - -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 face nodes */ - for(n1=0;n1vmapM[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;n1vmapP[id1] = k1*p_Np + mesh->Fmask[f1][n1]; - } - }else{ - /* treat as boundary for the moment */ - - for(n1=0;n1x[k1][mesh->Fmask[f1][n1]]; - y1 = mesh->y[k1][mesh->Fmask[f1][n1]]; - for(n2=0;n2x[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(d12vmapP[id1] = k2*p_Np + mesh->Fmask[f2][n2]; - } - } - } - } - } - } - -#if 0 - int n; - for(n=0;nK;++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;p2Npar[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;k1EToP[k1][f1]; - if(p2!=procid){ - for(n1=0;n1x[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;p2Npar[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;p2parNtotalout += 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;p2x[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(d12parmapOUT[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;p2EToP[k1][f1]==p2 && p2!=procid){ - for(n1=0;n1vmapP[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 face nodes */ - for(n1=0;n1vmapM[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;n1vmapP[id1] = k1*p_Np + mesh->Fmask[f1][n1]; - } - }else{ - /* treat as boundary for the moment */ - - for(n1=0;n1x[k1][mesh->Fmask[f1][n1]]; - y1 = mesh->y[k1][mesh->Fmask[f1][n1]]; - z1 = mesh->z[k1][mesh->Fmask[f1][n1]]; - - for(n2=0;n2x[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(d12vmapP[id1] = k2*p_Np + mesh->Fmask[f2][n2]; - break; - } - } - if(n2==p_Nfp){ - printf("LOST NODE !!!\n"); - } - } - } - } - } - -#if 0 - int n; - for(k1=0;k1K;++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;f1Nfaces;++f1){ - for(m=0;mx[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;p2Npar[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;k1EToP[k1][f1]; - if(p2!=procid){ - for(n1=0;n1x[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;p2Npar[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;p2parNtotalout += 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;p2x[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(d12parmapOUT[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;p2EToP[k1][f1]==p2 && p2!=procid){ - for(n1=0;n1vmapP[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(vb1vb2) - return 1; - else if(va1va2) - 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;kNpar = 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;nEToE[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;p2parK[p2] = BuildIntVector(mesh->Npar[p2]); - mesh->parF[p2] = BuildIntVector(mesh->Npar[p2]); - mesh->Npar[p2] = 0; - for(n=0;nparK[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(vc1vc2) - return 1; - else if(vb1vb2) - return 1; - else if(va1va2) - 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;kNpar = 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;nEToE[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;k1K;++k1){ - for(f1=0;f1Nfaces;++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;p2parK[p2] = BuildIntVector(mesh->Npar[p2]); - mesh->parF[p2] = BuildIntVector(mesh->Npar[p2]); - mesh->Npar[p2] = 0; - for(n=0;nparK[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;nf_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;nf_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;kK;++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;kK;++k){ - - GeometricFactors3d(mesh, k, - &drdx, &dsdx, &dtdx, - &drdy, &dsdy, &dtdy, - &drdz, &dsdz, &dtdz, &J); - - Normals3d(mesh, k, nxk, nyk, nzk, sJk); - - for(f=0;fNfaces;++f){ - - dt = min(dt, J/sJk[f]); - - for(m=0;mvmapM[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;kf_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;kf_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 - -#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;pGX!=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;pK*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;kK;++k){ - for(n=0;nx[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;nK*p_Np;++n) { - minEz = (minEz>Ez[n])?Ez[n]:minEz; - maxEz = (maxEz -#include - -texture t_LIFT; -texture t_DrDsDt; -texture t_Dr; -texture t_Ds; -texture t_Dt; -texture t_vgeo; -texture t_vgeo4; -texture t_Q; -texture t_partQ; -texture 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;mLIFT[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;nDr[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;kK;++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;kK;++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;fNfaces;++f){ - - dt = min(dt, J/sJk[f]); - - for(int m=0;mvmapM[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(nK; - 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>> (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 - -#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;pNpar[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;nparNtotalout;++n) - mesh->f_outQ[n] = f_Q[mesh->parmapOUT[n]]; - - /* do sends */ - int sk = 0, Nmess = 0; - for(p=0;pnprocs;++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;knprocs, sizeof(MPI_Status)); - MPI_Waitall(Nmess, mpi_in_requests, instatus); - free(instatus); - - for(k=0;knprocs, 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) { 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;nNv;++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;pK - nprocs*Klocal; - - - Kstart= 0; - for(p=0;pEToV = 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;nK;++n){ - fgets(buf, BUFSIZ, fp); - if(n>=Kstart && nEToV[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;vNverts;++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;nK;++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 -#include -#include -#include - -#include "mpi.h" - -#define max(a,b) ( (a>b)?a:b ) -#define min(a,b) ( (a */ - maxind = 0; - for(n=0;nr = BuildVector(p_Np); - mesh->s = BuildVector(p_Np); - mesh->t = BuildVector(p_Np); - for(n=0;nr[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;nDr[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;nLIFT[n][m] = p_LIFT[n][m]; - } - } - - mesh->Fmask = BuildIntMatrix(p_Nfaces, p_Nfp); - for(n=0;nFmask[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;kK;++k){ - for(n=0;nr[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;nGX[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 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 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; - -} - - -- cgit v1.3