| //------------------------------------------------------------------------------------------------------------------------------ |
| // Samuel Williams |
| // SWWilliams@lbl.gov |
| // Lawrence Berkeley National Lab |
| //------------------------------------------------------------------------------------------------------------------------------ |
| #include <stdint.h> |
| #include "timer.h" |
| #include "mg.h" |
| //------------------------------------------------------------------------------------------------------------------------------ |
| // piecewise constant interpolation... |
| // |
| // +-------+ +---+---+ |
| // | | | x | x | |
| // | x | -> +---+---+ |
| // | | | x | x | |
| // +-------+ +---+---+ |
| // |
| //------------------------------------------------------------------------------------------------------------------------------ |
| void interpolation_constant(domain_type * domain, int level_f, double prescale_f, int id_f, int id_c){ // id_f = prescale*id_f + I_{2h}^{h}(id_c) |
| int level_c = level_f+1; |
| uint64_t _timeStart = CycleTime(); |
| |
| int CollaborativeThreadingBoxSize = 100000; // i.e. never |
| #ifdef __COLLABORATIVE_THREADING |
| CollaborativeThreadingBoxSize = 1 << __COLLABORATIVE_THREADING; |
| #endif |
| int omp_across_boxes = (domain->subdomains[0].levels[level_f].dim.i < CollaborativeThreadingBoxSize); |
| int omp_within_a_box = (domain->subdomains[0].levels[level_f].dim.i >= CollaborativeThreadingBoxSize); |
| int box; |
| |
| #pragma omp parallel for private(box) if(omp_across_boxes) |
| for(box=0;box<domain->subdomains_per_rank;box++){ |
| int i,j,k; |
| int ghosts_c = domain->subdomains[box].levels[level_c].ghosts; |
| int pencil_c = domain->subdomains[box].levels[level_c].pencil; |
| int plane_c = domain->subdomains[box].levels[level_c].plane; |
| |
| int ghosts_f = domain->subdomains[box].levels[level_f].ghosts; |
| int pencil_f = domain->subdomains[box].levels[level_f].pencil; |
| int plane_f = domain->subdomains[box].levels[level_f].plane; |
| int dim_i_f = domain->subdomains[box].levels[level_f].dim.i; |
| int dim_j_f = domain->subdomains[box].levels[level_f].dim.j; |
| int dim_k_f = domain->subdomains[box].levels[level_f].dim.k; |
| |
| double * __restrict__ grid_f = domain->subdomains[box].levels[level_f].grids[id_f] + ghosts_f*(1+pencil_f+plane_f); |
| double * __restrict__ grid_c = domain->subdomains[box].levels[level_c].grids[id_c] + ghosts_c*(1+pencil_c+plane_c); |
| |
| #pragma omp parallel for private(k,j,i) if(omp_within_a_box) collapse(2) |
| for(k=0;k<dim_k_f;k++){ |
| for(j=0;j<dim_j_f;j++){ |
| for(i=0;i<dim_i_f;i++){ |
| int ijk_f = (i ) + (j )*pencil_f + (k )*plane_f; |
| int ijk_c = (i>>1) + (j>>1)*pencil_c + (k>>1)*plane_c; |
| grid_f[ijk_f] = prescale_f*grid_f[ijk_f] + grid_c[ijk_c]; |
| }}} |
| } |
| domain->cycles.interpolation[level_f] += (uint64_t)(CycleTime()-_timeStart); |
| } |
| |
| //------------------------------------------------------------------------------------------------------------------------------ |
| // piecewise linear interpolation... |
| // |
| //------------------------------------------------------------------------------------------------------------------------------ |
| void interpolation_linear(domain_type * domain, int level_f, double prescale_f, int id_f, int id_c){ // id_f = prescale*id_f + I_{2h}^{h}(id_c) |
| int level_c = level_f+1; |
| exchange_boundary(domain,level_c,id_c,1,1,1); // linear needs corners/edges in the coarse grid. |
| |
| uint64_t _timeStart = CycleTime(); |
| |
| int CollaborativeThreadingBoxSize = 100000; // i.e. never |
| #ifdef __COLLABORATIVE_THREADING |
| CollaborativeThreadingBoxSize = 1 << __COLLABORATIVE_THREADING; |
| #endif |
| int omp_across_boxes = (domain->subdomains[0].levels[level_f].dim.i < CollaborativeThreadingBoxSize); |
| int omp_within_a_box = (domain->subdomains[0].levels[level_f].dim.i >= CollaborativeThreadingBoxSize); |
| int box; |
| |
| #pragma omp parallel for private(box) if(omp_across_boxes) |
| for(box=0;box<domain->subdomains_per_rank;box++){ |
| int i,j,k; |
| int ghosts_c = domain->subdomains[box].levels[level_c].ghosts; |
| int pencil_c = domain->subdomains[box].levels[level_c].pencil; |
| int plane_c = domain->subdomains[box].levels[level_c].plane; |
| int dim_i_c = domain->subdomains[box].levels[level_c].dim.i; |
| int dim_j_c = domain->subdomains[box].levels[level_c].dim.j; |
| int dim_k_c = domain->subdomains[box].levels[level_c].dim.k; |
| |
| int ghosts_f = domain->subdomains[box].levels[level_f].ghosts; |
| int pencil_f = domain->subdomains[box].levels[level_f].pencil; |
| int plane_f = domain->subdomains[box].levels[level_f].plane; |
| int dim_i_f = domain->subdomains[box].levels[level_f].dim.i; |
| int dim_j_f = domain->subdomains[box].levels[level_f].dim.j; |
| int dim_k_f = domain->subdomains[box].levels[level_f].dim.k; |
| |
| double * __restrict__ grid_f = domain->subdomains[box].levels[level_f].grids[ id_f] + ghosts_f*(1 + pencil_f + plane_f); // [0] is first non-ghost zone element |
| double * __restrict__ grid_c = domain->subdomains[box].levels[level_c].grids[ id_c] + ghosts_c*(1 + pencil_c + plane_c); |
| |
| // FIX what about dirichlet boundary conditions ??? |
| #pragma omp parallel for private(k,j,i) if(omp_within_a_box) collapse(2) |
| for(k=0;k<dim_k_f;k++){ |
| for(j=0;j<dim_j_f;j++){ |
| for(i=0;i<dim_i_f;i++){ |
| int ijk_f = (i ) + (j )*pencil_f + (k )*plane_f; |
| int ijk_c = (i>>1) + (j>>1)*pencil_c + (k>>1)*plane_c; |
| |
| // ----------------------------------------------------------------------------------------------------------------------- |
| // Piecewise Quadratic Interpolation |
| // ----------------------------------------------------------------------------------------------------------------------- |
| // define parabola f(x) = ax^2 + bx + c through three coarse grid cells in i-direction... x=-1, x=0, and x=1 |
| // interpolate to (-0.25,j,k) and (+0.25,j,k) |
| // combine into 3 dimensions |
| // |
| // +-------+-------+-------+ |
| // | o | o | o | |
| // +-------+-------+-------+ |
| // . |
| // . |
| // interpolation |
| // . |
| // . |
| // +---+---+---+---+---+---+ |
| // | | | x | y | | | |
| // +---+---+---+---+---+---+ |
| // |
| #warning using 27pt stencil for piecewise-quadratic interpolation |
| double xm= 0.156250,x0=0.937500,xp=-0.093750; |
| double ym= 0.156250,y0=0.937500,yp=-0.093750; |
| double zm= 0.156250,z0=0.937500,zp=-0.093750; |
| if(i&0x1){xm=-0.093750;x0=0.937500;xp= 0.156250;} |
| if(j&0x1){ym=-0.093750;y0=0.937500;yp= 0.156250;} |
| if(k&0x1){zm=-0.093750;z0=0.937500;zp= 0.156250;} |
| grid_f[ijk_f] = |
| prescale_f*grid_f[ijk_f ] + |
| |
| zm*ym*xm*grid_c[ijk_c-1-pencil_c-plane_c] + |
| zm*ym*x0*grid_c[ijk_c -pencil_c-plane_c] + |
| zm*ym*xp*grid_c[ijk_c+1-pencil_c-plane_c] + |
| zm*y0*xm*grid_c[ijk_c-1 -plane_c] + |
| zm*y0*x0*grid_c[ijk_c -plane_c] + |
| zm*y0*xp*grid_c[ijk_c+1 -plane_c] + |
| zm*yp*xm*grid_c[ijk_c-1+pencil_c-plane_c] + |
| zm*yp*x0*grid_c[ijk_c +pencil_c-plane_c] + |
| zm*yp*xp*grid_c[ijk_c+1+pencil_c-plane_c] + |
| |
| z0*ym*xm*grid_c[ijk_c-1-pencil_c ] + |
| z0*ym*x0*grid_c[ijk_c -pencil_c ] + |
| z0*ym*xp*grid_c[ijk_c+1-pencil_c ] + |
| z0*y0*xm*grid_c[ijk_c-1 ] + |
| z0*y0*x0*grid_c[ijk_c ] + |
| z0*y0*xp*grid_c[ijk_c+1 ] + |
| z0*yp*xm*grid_c[ijk_c-1+pencil_c ] + |
| z0*yp*x0*grid_c[ijk_c +pencil_c ] + |
| z0*yp*xp*grid_c[ijk_c+1+pencil_c ] + |
| |
| zp*ym*xm*grid_c[ijk_c-1-pencil_c+plane_c] + |
| zp*ym*x0*grid_c[ijk_c -pencil_c+plane_c] + |
| zp*ym*xp*grid_c[ijk_c+1-pencil_c+plane_c] + |
| zp*y0*xm*grid_c[ijk_c-1 +plane_c] + |
| zp*y0*x0*grid_c[ijk_c +plane_c] + |
| zp*y0*xp*grid_c[ijk_c+1 +plane_c] + |
| zp*yp*xm*grid_c[ijk_c-1+pencil_c+plane_c] + |
| zp*yp*x0*grid_c[ijk_c +pencil_c+plane_c] + |
| zp*yp*xp*grid_c[ijk_c+1+pencil_c+plane_c] ; |
| /* |
| // ----------------------------------------------------------------------------------------------------------------------- |
| // Piecewise-Linear Interpolation |
| // ----------------------------------------------------------------------------------------------------------------------- |
| // define line f(x) = bx + c through two coarse grid cells in i-direction... x=+/-1 and x=0 |
| // interpolate to either (-0.25) or (+0.25) |
| // combine into 3 dimensions |
| // |
| // +-------+-------+ |
| // | o | o | |
| // +-------+-------+ |
| // . |
| // . |
| // interpolation |
| // . |
| // . |
| // +---+---+---+---+ |
| // | | x | y | | |
| // +---+---+---+---+ |
| // |
| #warning using 8pt stencil for piecewise-linear interpolation |
| int delta_i= -1;if(i&0x1)delta_i= 1; // i.e. even points look backwards while odd points look forward |
| int delta_j=-pencil_c;if(j&0x1)delta_j=pencil_c; |
| int delta_k= -plane_c;if(k&0x1)delta_k= plane_c; |
| grid_f[ijk_f] = |
| prescale_f*grid_f[ijk_f ] + |
| 0.421875*grid_c[ijk_c ] + |
| 0.140625*grid_c[ijk_c +delta_k] + |
| 0.140625*grid_c[ijk_c +delta_j ] + |
| 0.046875*grid_c[ijk_c +delta_j+delta_k] + |
| 0.140625*grid_c[ijk_c+delta_i ] + |
| 0.046875*grid_c[ijk_c+delta_i +delta_k] + |
| 0.046875*grid_c[ijk_c+delta_i+delta_j ] + |
| 0.015625*grid_c[ijk_c+delta_i+delta_j+delta_k] ; |
| */ |
| /* |
| // ----------------------------------------------------------------------------------------------------------------------- |
| // Piecewise-Linear Interpolation |
| // ----------------------------------------------------------------------------------------------------------------------- |
| #warning using 7pt stencil for piecewise-linear interpolation... doesn't given 2nd order FMG?? |
| double coefi = -0.125;if(i&0x1)coefi = 0.125; |
| double coefj = -0.125;if(j&0x1)coefj = 0.125; |
| double coefk = -0.125;if(k&0x1)coefk = 0.125; |
| grid_f[ijk_f] = prescale_f*grid_f[ijk_f ] + |
| grid_c[ijk_c ] + |
| coefi*(grid_c[ijk_c+ 1]-grid_c[ijk_c- 1]) + |
| coefj*(grid_c[ijk_c+pencil_c]-grid_c[ijk_c-pencil_c]) + |
| coefk*(grid_c[ijk_c+ plane_c]-grid_c[ijk_c- plane_c]); |
| */ |
| }}} |
| } |
| domain->cycles.interpolation[level_f] += (uint64_t)(CycleTime()-_timeStart); |
| } |
| //------------------------------------------------------------------------------------------------------------------------------ |
