diff options
Diffstat (limited to 'benchmarks/CUDA/DG/src/InitCPU3d.c')
| -rw-r--r-- | benchmarks/CUDA/DG/src/InitCPU3d.c | 152 |
1 files changed, 152 insertions, 0 deletions
diff --git a/benchmarks/CUDA/DG/src/InitCPU3d.c b/benchmarks/CUDA/DG/src/InitCPU3d.c new file mode 100644 index 0000000..2fe907f --- /dev/null +++ b/benchmarks/CUDA/DG/src/InitCPU3d.c @@ -0,0 +1,152 @@ +#include "mpi.h" +#include "fem.h" + +double InitCPU3d(Mesh *mesh, int Nfields){ + + printf("Np = %d, BSIZE = %d\n", p_Np, BSIZE); + + /* Q */ + int sz = mesh->K*(p_Np)*Nfields*sizeof(float); /* TW BLOCK */ + + mesh->f_Q = (float*) calloc(mesh->K*p_Np*Nfields, sizeof(float)); + mesh->f_rhsQ = (float*) calloc(mesh->K*p_Np*Nfields, sizeof(float)); + mesh->f_resQ = (float*) calloc(mesh->K*p_Np*Nfields, sizeof(float)); + + /* float LIFT */ + sz = p_Np*(p_Nfp)*(p_Nfaces)*sizeof(float); + mesh->f_LIFT = (float*) malloc(sz); + int sk = 0, n, m, f, k; + + for(n=0;n<p_Np;++n){ + for(m=0;m<p_Nfp*p_Nfaces;++m){ + mesh->f_LIFT[sk++] = mesh->LIFT[n][m]; + } + } + + /* float Dr & Ds */ + sz = p_Np*p_Np*sizeof(float); + mesh->f_Dr = (float*) malloc(sz); + mesh->f_Ds = (float*) malloc(sz); + mesh->f_Dt = (float*) malloc(sz); + + sk = 0; + for(n=0;n<p_Np;++n){ + for(m=0;m<p_Np;++m){ + mesh->f_Dr[sk] = mesh->Dr[n][m]; + mesh->f_Ds[sk] = mesh->Ds[n][m]; + mesh->f_Dt[sk] = mesh->Dt[n][m]; + ++sk; + } + } + + /* vgeo */ + double drdx, dsdx, dtdx; + double drdy, dsdy, dtdy; + double drdz, dsdz, dtdz, J; + mesh->vgeo = (float*) calloc(12*mesh->K, sizeof(float)); + + for(k=0;k<mesh->K;++k){ + GeometricFactors3d(mesh, k, + &drdx, &dsdx, &dtdx, + &drdy, &dsdy, &dtdy, + &drdz, &dsdz, &dtdz, &J); + + mesh->vgeo[k*12+0] = drdx; mesh->vgeo[k*12+1] = drdy; mesh->vgeo[k*12+2] = drdz; + mesh->vgeo[k*12+4] = dsdx; mesh->vgeo[k*12+5] = dsdy; mesh->vgeo[k*12+6] = dsdz; + mesh->vgeo[k*12+8] = dtdx; mesh->vgeo[k*12+9] = dtdy; mesh->vgeo[k*12+10] = dtdz; + } + + /* surfinfo (vmapM, vmapP, Fscale, Bscale, nx, ny, nz, 0) */ + sz = mesh->K*p_Nfp*p_Nfaces*7*sizeof(float); + + mesh->surfinfo = (float*) malloc(sz); + + /* local-local info */ + sk = 0; + int skP = -1; + double *nxk = BuildVector(mesh->Nfaces); + double *nyk = BuildVector(mesh->Nfaces); + double *nzk = BuildVector(mesh->Nfaces); + double *sJk = BuildVector(mesh->Nfaces); + + double dt = 1e6; + + sk = 0; + for(k=0;k<mesh->K;++k){ + + GeometricFactors3d(mesh, k, + &drdx, &dsdx, &dtdx, + &drdy, &dsdy, &dtdy, + &drdz, &dsdz, &dtdz, &J); + + Normals3d(mesh, k, nxk, nyk, nzk, sJk); + + for(f=0;f<mesh->Nfaces;++f){ + + dt = min(dt, J/sJk[f]); + + for(m=0;m<p_Nfp;++m){ + int id = m + f*p_Nfp + p_Nfp*p_Nfaces*k; + int idM = mesh->vmapM[id]; + int idP = mesh->vmapP[id]; + int nM = idM%p_Np; + int nP = idP%p_Np; + int kM = (idM-nM)/p_Np; + int kP = (idP-nP)/p_Np; + idM = Nfields*(nM + p_Np*kM); + idP = Nfields*(nP + p_Np*kP); + + /* stub resolve some other way */ + if(mesh->vmapP[id]<0){ + idP = mesh->vmapP[id]; /* -ve numbers */ + } + + mesh->surfinfo[sk++] = idM; + mesh->surfinfo[sk++] = idP; + mesh->surfinfo[sk++] = sJk[f]/(2.*J); + mesh->surfinfo[sk++] = (idM==idP)?-1.:1.; + mesh->surfinfo[sk++] = nxk[f]; + mesh->surfinfo[sk++] = nyk[f]; + mesh->surfinfo[sk++] = nzk[f]; + } + } + } +} + +void cpu_set_data3d(Mesh *mesh, double *Hx, double *Hy, double *Hz, + double *Ex, double *Ey, double *Ez){ + + const int K = mesh->K; + int k, n, sk=0; + + for(k=0;k<K;++k){ + for(n=0;n<p_Np;++n){ + mesh->f_Q[sk++] = Hx[n+k*p_Np]; + mesh->f_Q[sk++] = Hy[n+k*p_Np]; + mesh->f_Q[sk++] = Hz[n+k*p_Np]; + mesh->f_Q[sk++] = Ex[n+k*p_Np]; + mesh->f_Q[sk++] = Ey[n+k*p_Np]; + mesh->f_Q[sk++] = Ez[n+k*p_Np]; + } + } + +} + + +void cpu_get_data3d(Mesh *mesh, double *Hx, double *Hy, double *Hz, + double *Ex, double *Ey, double *Ez){ + const int K = mesh->K; + int k, n, sk=0; + + for(k=0;k<K;++k){ + for(n=0;n<p_Np;++n){ + Hx[n+k*p_Np] = mesh->f_Q[sk++]; + Hy[n+k*p_Np] = mesh->f_Q[sk++]; + Hz[n+k*p_Np] = mesh->f_Q[sk++]; + Ex[n+k*p_Np] = mesh->f_Q[sk++]; + Ey[n+k*p_Np] = mesh->f_Q[sk++]; + Ez[n+k*p_Np] = mesh->f_Q[sk++]; + } + } + +} |