From 69f2911e04ffb1b19eef1fafb8c040af271f656e Mon Sep 17 00:00:00 2001 From: Tor Aamodt Date: Thu, 15 Jul 2010 18:09:46 -0800 Subject: creating branch for adding support for CUDA 3.x and Fermi [git-p4: depot-paths = "//depot/gpgpu_sim_research/fermi/distribution/": change = 6829] --- 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 insertions(+) create mode 100644 benchmarks/CUDA/DG/src/BuildMaps2d.c create mode 100644 benchmarks/CUDA/DG/src/BuildMaps3d.c create mode 100644 benchmarks/CUDA/DG/src/FacePair2d.c create mode 100644 benchmarks/CUDA/DG/src/FacePair3d.c create mode 100644 benchmarks/CUDA/DG/src/InitCPU3d.c create mode 100644 benchmarks/CUDA/DG/src/LoadBalance3d.c create mode 100644 benchmarks/CUDA/DG/src/MaxwellsDriver3d.c create mode 100644 benchmarks/CUDA/DG/src/MaxwellsKernel3d.cu create mode 100644 benchmarks/CUDA/DG/src/MaxwellsMPI3d.c create mode 100644 benchmarks/CUDA/DG/src/MaxwellsRHS3d.c create mode 100644 benchmarks/CUDA/DG/src/MaxwellsRun3d.c create mode 100644 benchmarks/CUDA/DG/src/Mesh3d.c create mode 100644 benchmarks/CUDA/DG/src/ParallelPairs.c create mode 100644 benchmarks/CUDA/DG/src/StartUp3d.c create mode 100644 benchmarks/CUDA/DG/src/Utils.c create mode 100644 benchmarks/CUDA/DG/src/buildsource.m create mode 100644 benchmarks/CUDA/DG/src/gpuDG.vpj create 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 new file mode 100644 index 0000000..7d0d57c --- /dev/null +++ b/benchmarks/CUDA/DG/src/BuildMaps2d.c @@ -0,0 +1,231 @@ +#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 new file mode 100644 index 0000000..89f7a44 --- /dev/null +++ b/benchmarks/CUDA/DG/src/BuildMaps3d.c @@ -0,0 +1,279 @@ +#include "mpi.h" +#include "fem.h" + +void BuildMaps3d(Mesh *mesh){ + + printf("Hello %d\n", 1002); + + int nprocs = mesh->nprocs; + int procid = mesh->procid; + + int K = mesh->K; + int Nfaces = mesh->Nfaces; + + mesh->vmapM = BuildIntVector(p_Nfp*p_Nfaces*K); + mesh->vmapP = BuildIntVector(p_Nfp*p_Nfaces*K); + + int m; + int k1,f1,p1,n1,id1, k2,f2,p2,n2,id2; + + double x1, y1,z1, x2, y2, z2, d12; + + double *nxk = BuildVector(Nfaces); + double *nyk = BuildVector(Nfaces); + double *nzk = BuildVector(Nfaces); + double *sJk = BuildVector(Nfaces); + + printf("Hello %d\n", 1001); + + /* first build local */ + for(k1=0;k1 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 new file mode 100644 index 0000000..92e47f5 --- /dev/null +++ b/benchmarks/CUDA/DG/src/FacePair2d.c @@ -0,0 +1,154 @@ +#include "mpi.h" +#include "fem.h" + +typedef struct foob { + int p1, k1, f1, p2, k2, f2, va, vb, g; +}face; + +int compare_pairs(const void *obj1, const void *obj2){ + + face *e1 = (face*) obj1; + face *e2 = (face*) obj2; + + int a1 = e1->va, b1 = e1->vb; + int a2 = e2->va, b2 = e2->vb; + + int va1, vb1, va2, vb2; + + va1 = min(a1, b1); + vb1 = max(a1, b1); + + va2 = min(a2, b2); + vb2 = max(a2, b2); + + if(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 new file mode 100644 index 0000000..41d1fad --- /dev/null +++ b/benchmarks/CUDA/DG/src/FacePair3d.c @@ -0,0 +1,201 @@ +#include "mpi.h" +#include "fem.h" + +typedef struct foob { + int p1, k1, f1, p2, k2, f2; + int va, vb, vc, g; +}face3d; + +int compare_pairs3d(const void *obj1, const void *obj2){ + + face3d *e1 = (face3d*) obj1; + face3d *e2 = (face3d*) obj2; + + int a1 = e1->va, b1 = e1->vb, c1 = e1->vc; + int a2 = e2->va, b2 = e2->vb, c2 = e2->vc; + + int va1, vb1, vc1, va2, vb2, vc2; + + va1 = min(a1, min(b1, c1)); + vc1 = max(a1, max(b1, c1)); + + if(va1!=a1 && vc1!=a1) vb1=a1; + if(va1!=b1 && vc1!=b1) vb1=b1; + if(va1!=c1 && vc1!=c1) vb1=c1; + + va2 = min(a2, min(b2, c2)); + vc2 = max(a2, max(b2, c2)); + + if(va2!=a2 && vc2!=a2) vb2=a2; + if(va2!=b2 && vc2!=b2) vb2=b2; + if(va2!=c2 && vc2!=c2) vb2=c2; + + if(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 new file mode 100644 index 0000000..2fe907f --- /dev/null +++ b/benchmarks/CUDA/DG/src/InitCPU3d.c @@ -0,0 +1,152 @@ +#include "mpi.h" +#include "fem.h" + +double InitCPU3d(Mesh *mesh, int Nfields){ + + printf("Np = %d, BSIZE = %d\n", p_Np, BSIZE); + + /* Q */ + int sz = mesh->K*(p_Np)*Nfields*sizeof(float); /* TW BLOCK */ + + mesh->f_Q = (float*) calloc(mesh->K*p_Np*Nfields, sizeof(float)); + mesh->f_rhsQ = (float*) calloc(mesh->K*p_Np*Nfields, sizeof(float)); + mesh->f_resQ = (float*) calloc(mesh->K*p_Np*Nfields, sizeof(float)); + + /* float LIFT */ + sz = p_Np*(p_Nfp)*(p_Nfaces)*sizeof(float); + mesh->f_LIFT = (float*) malloc(sz); + int sk = 0, n, m, f, k; + + for(n=0;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 new file mode 100644 index 0000000..586cd8c --- /dev/null +++ b/benchmarks/CUDA/DG/src/LoadBalance3d.c @@ -0,0 +1,243 @@ +#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 new file mode 100644 index 0000000..6bdf46f --- /dev/null +++ b/benchmarks/CUDA/DG/src/MaxwellsMPI3d.c @@ -0,0 +1,71 @@ + +#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 new file mode 100644 index 0000000..13ed5bf --- /dev/null +++ b/benchmarks/CUDA/DG/src/MaxwellsRHS3d.c @@ -0,0 +1,220 @@ +#include "mpi.h" +#include "fem.h" + +void MaxwellsRHS3d(Mesh *mesh, float frka, float frkb, float fdt){ + + /* registers and temporary */ + register unsigned int k, n; + + /* mesh parameters */ + const int K = mesh->K; + + float *vgeo = mesh->vgeo; + float *surfinfo = mesh->surfinfo; + float *f_Dr = mesh->f_Dr; + float *f_Ds = mesh->f_Ds; + float *f_Dt = mesh->f_Dt; + float *f_LIFT = mesh->f_LIFT; + + float *f_Q = mesh->f_Q; + float *f_rhsQ = mesh->f_rhsQ; + float *f_resQ = mesh->f_resQ; + + float *f_inQ = mesh->f_inQ; + float *f_outQ = mesh->f_outQ; + + int p; + + /* mpi request buffer */ + MPI_Request *mpi_out_requests = (MPI_Request*) calloc(mesh->nprocs, sizeof(MPI_Request)); + MPI_Request *mpi_in_requests = (MPI_Request*) calloc(mesh->nprocs, sizeof(MPI_Request)); + + /* buffer outgoing node data */ + for(n=0;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 new file mode 100644 index 0000000..d4ee9f9 --- /dev/null +++ b/benchmarks/CUDA/DG/src/MaxwellsRun3d.c @@ -0,0 +1,60 @@ +#include "mpi.h" +#include "fem.h" + +void MaxwellsRun3d(Mesh *mesh, double FinalTime, double dt){ + + double time = 0; + int INTRK, tstep=0; + + double mpitime0 = MPI_Wtime(); + MPI_Barrier(MPI_COMM_WORLD); + + dt = 0.001; + + /* outer time step loop */ + while (time FinalTime) { dt = FinalTime-time; } + + for (INTRK=1; INTRK<=5; ++INTRK) { + + /* compute rhs of TM-mode MaxwellsGPU's equations */ + const float fdt = dt; + const float fa = (float)mesh->rk4a[INTRK-1]; + const float fb = (float)mesh->rk4b[INTRK-1]; + +#ifdef CUDA + MaxwellsKernel3d(mesh, fa, fb, fdt); +#else + MaxwellsRHS3d(mesh, fa, fb, fdt); +#endif + } + + time += dt; /* increment current time */ + tstep++; /* increment timestep */ + } + +#ifdef CUDA + cudaThreadSynchronize(); +#endif + + double flopsV = p_Np*p_Np*36 + p_Np*66; /* V3 */ + double flopsS = p_Nfp*p_Nfaces*(15 + 10 + 36) + p_Np*(p_Nfaces*p_Nfp*12 + 6); + double flopsR = p_Np*p_Nfields*4; + + int Kloc = mesh->K; + + double mpitime1 = MPI_Wtime(); + + double time_total = mpitime1-mpitime0; + + MPI_Barrier(MPI_COMM_WORLD); + printf("%d %d %lg %lg proc, N, time taken, GFLOPS/s (GPU)\n", + mesh->procid, + p_N, + time_total, + 5*(tstep)*( (flopsV+flopsS+flopsR)*((double)Kloc/(1.e9*time_total)))); + + +} diff --git a/benchmarks/CUDA/DG/src/Mesh3d.c b/benchmarks/CUDA/DG/src/Mesh3d.c new file mode 100644 index 0000000..af76acb --- /dev/null +++ b/benchmarks/CUDA/DG/src/Mesh3d.c @@ -0,0 +1,195 @@ +#include "mpi.h" +#include "fem.h" + +Mesh *ReadMesh3d(char *filename){ + + int n; + + Mesh *mesh = (Mesh*) calloc(1, sizeof(Mesh)); + + char buf[BUFSIZ]; + + FILE *fp = fopen(filename, "r"); + + /* assume modified Gambit neutral format */ + for(n=0;n<6;++n) + fgets(buf, BUFSIZ, fp); + + fgets(buf, BUFSIZ, fp); + sscanf(buf, "%d %d \n", &(mesh->Nv), &(mesh->K)); + mesh->Nverts = 4; /* assume tets */ + mesh->Nedges = 6; /* assume tets */ + mesh->Nfaces = 4; /* assume tets */ + + fgets(buf, BUFSIZ, fp); + fgets(buf, BUFSIZ, fp); + + /* read vertex coordinates */ + double *VX = BuildVector(mesh->Nv); + double *VY = BuildVector(mesh->Nv); + double *VZ = BuildVector(mesh->Nv); + for(n=0;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 new file mode 100644 index 0000000..7474c0b --- /dev/null +++ b/benchmarks/CUDA/DG/src/ParallelPairs.c @@ -0,0 +1,124 @@ +#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 new file mode 100644 index 0000000..e0364cf --- /dev/null +++ b/benchmarks/CUDA/DG/src/Utils.c @@ -0,0 +1,155 @@ +#include "fem.h" + +/* some very basic memory allocation routines */ + +/* row major storage for a 2D matrix array */ +double **BuildMatrix(int Nrows, int Ncols){ + int n; + double **A = (double**) calloc(Nrows, sizeof(double*)); + + A[0] = (double*) calloc(Nrows*Ncols, sizeof(double)); + + for(n=1;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 new file mode 100644 index 0000000..c93d363 --- /dev/null +++ b/benchmarks/CUDA/DG/src/buildsource.m @@ -0,0 +1,106 @@ + +function buildsource(in_N) + +Globals2D; + + +NODETOL = 1e-8; + +N = in_N; + +Np = (N+1)*(N+2)/2; Nfp = N+1; Nfaces = 3; + +[x,y] = Nodes2D(N); +[r,s] = xytors(x,y); + +% find all the nodes that lie on each edge +fmask1 = find( abs(s+1) < NODETOL)'; +fmask2 = find( abs(r+s) < NODETOL)'; +fmask3 = find( abs(r+1) < NODETOL)'; +Fmask = [fmask1;fmask2;fmask3]'; + +% vandermonde matrix +V = Vandermonde2D(N, r, s); + +% derivative matrices +[Dr, Ds] = Dmatrices2D(N, r, s, V); + +% lift matrix +LIFT = Lift2D(); + +% correct to 0-index +Fmask = Fmask-1; + + +fid = fopen(sprintf('dataN%02d.h', N), 'w'); + +fprintf(fid, sprintf('#ifndef DATAN%02d \n\n', N)); +fprintf(fid, sprintf('#define DATAN%02d 1\n\n', N)); +fprintf(fid, 'double p_r[%d] = { ', Np); +fprintf(fid, '%17.15g ', r(1)); +for n=2:Np + fprintf(fid, ', %17.15g ', r(n)); +end +fprintf(fid, '};\n') + +fprintf(fid, 'double p_s[%d] = {', Np); +fprintf(fid, '%17.15g ', s(1)); +for n=2:Np + fprintf(fid, ', %17.15g ', s(n)); +end +fprintf(fid, '};\n') + +fprintf(fid, 'double p_Dr[%d][%d] = {', Np, Np); +for m=1:Np +fprintf(fid, '{%17.15g ', Dr(m,1)); +for n=2:Np + fprintf(fid, ', %17.15g ', Dr(m,n)); +end +if(m + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/benchmarks/CUDA/DG/src/tictoc.cu b/benchmarks/CUDA/DG/src/tictoc.cu new file mode 100644 index 0000000..f1335c3 --- /dev/null +++ b/benchmarks/CUDA/DG/src/tictoc.cu @@ -0,0 +1,23 @@ +cudaEvent_t cstart, cstop; + +void cudatic(){ + + cudaEventCreate(&cstart); + cudaEventCreate(&cstop); + + cudaEventRecord(cstart, 0); +} + +float cudatoc(){ + + cudaEventRecord(cstop, 0); + cudaEventSynchronize(cstop); + float elapsedTime; + cudaEventElapsedTime(&elapsedTime, cstart, cstop); + + /* return elapsed time in seconds */ + return elapsedTime/1000.0; + +} + + -- cgit v1.3