in MultiSource/Benchmarks/DOE-ProxyApps-C++/CLAMR/mesh.cpp [4220:6277]
void Mesh::calc_neighbors_local(void)
{
struct timeval tstart_cpu;
cpu_timer_start(&tstart_cpu);
if (do_rezone) {
int flags = INDEX_ARRAY_MEMORY;
#if defined (HAVE_J7)
if (parallel) flags |= LOAD_BALANCE_MEMORY;
#endif
#ifdef _OPENMP
#pragma omp master
{
#endif
cpu_counters[MESH_COUNTER_CALC_NEIGH]++;
if (mesh_memory.get_memory_size(nlft) < ncells){
if (nlft != NULL) nlft = (int *)mesh_memory.memory_delete(nlft);
if (nrht != NULL) nrht = (int *)mesh_memory.memory_delete(nrht);
if (nbot != NULL) nbot = (int *)mesh_memory.memory_delete(nbot);
if (ntop != NULL) ntop = (int *)mesh_memory.memory_delete(ntop);
nlft = (int *)mesh_memory.memory_malloc(ncells, sizeof(int), "nlft", flags);
nrht = (int *)mesh_memory.memory_malloc(ncells, sizeof(int), "nrht", flags);
nbot = (int *)mesh_memory.memory_malloc(ncells, sizeof(int), "nbot", flags);
ntop = (int *)mesh_memory.memory_malloc(ncells, sizeof(int), "ntop", flags);
}
#ifdef _OPENMP
}
#pragma omp barrier
#endif
int lowerBound, upperBound;
set_bounds(ncells);
get_bounds(lowerBound, upperBound);
for (int ic = lowerBound; ic < upperBound; ic++){
nlft[ic] = -98;
nrht[ic] = -98;
nbot[ic] = -98;
ntop[ic] = -98;
}
if (calc_neighbor_type == HASH_TABLE) {
struct timeval tstart_lev2;
if (TIMING_LEVEL >= 2) cpu_timer_start(&tstart_lev2);
ncells_ghost = ncells;
// Find maximum i column and j row for this processor
static int jmintile, imintile, jmaxtile, imaxtile;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
jmintile = (jmax+1)*IPOW2(levmx);
imintile = (imax+1)*IPOW2(levmx);
jmaxtile = 0;
imaxtile = 0;
#ifdef _OPENMP
}
#pragma omp barrier
#endif
int my_jmintile = jmintile;
int my_imintile = imintile;
int my_jmaxtile = 0;
int my_imaxtile = 0;
#ifdef _OPENMP
#pragma omp for
#endif
for(uint ic=0; ic<ncells; ic++){
int lev = level[ic];
// if (lev < 0 || lev > levmx) printf("DEBUG -- cell %d lev %d\n",ic,level[ic]);
if ( j[ic] *IPOW2(levmx-lev) < my_jmintile) my_jmintile = j[ic] *IPOW2(levmx-lev) ;
if ((j[ic]+1)*IPOW2(levmx-lev)-1 > my_jmaxtile) my_jmaxtile = (j[ic]+1)*IPOW2(levmx-lev)-1;
if ( i[ic] *IPOW2(levmx-lev) < my_imintile) my_imintile = i[ic] *IPOW2(levmx-lev) ;
if ((i[ic]+1)*IPOW2(levmx-lev)-1 > my_imaxtile) my_imaxtile = (i[ic]+1)*IPOW2(levmx-lev)-1;
}
#ifdef _OPENMP
#pragma omp critical
{
#endif
if (my_jmintile < jmintile) jmintile = my_jmintile;
if (my_imintile < imintile) imintile = my_imintile;
if (my_jmaxtile > jmaxtile) jmaxtile = my_jmaxtile;
if (my_imaxtile > imaxtile) imaxtile = my_imaxtile;
#ifdef _OPENMP
} // end critical region
#pragma omp barrier
#endif
//if (DEBUG) fprintf(fp,"%d: Tile Sizes are imin %d imax %d jmin %d jmax %d\n",mype,imintile,imaxtile,jmintile,jmaxtile);
// Expand size by 2*coarse_cells for ghost cells
int jminsize = max(jmintile-2*IPOW2(levmx),0);
int jmaxsize = min(jmaxtile+2*IPOW2(levmx),(jmax+1)*IPOW2(levmx));
int iminsize = max(imintile-2*IPOW2(levmx),0);
int imaxsize = min(imaxtile+2*IPOW2(levmx),(imax+1)*IPOW2(levmx));
//if (DEBUG) fprintf(fp,"%d: Sizes are imin %d imax %d jmin %d jmax %d\n",mype,iminsize,imaxsize,jminsize,jmaxsize);
//fprintf(fp,"DEBUG -- ncells %lu\n",ncells);
static int *hash;
#ifdef _OPENMP
hash = compact_hash_init_openmp(ncells, imaxsize-iminsize, jmaxsize-jminsize, 0);
#else
hash = compact_hash_init(ncells, imaxsize-iminsize, jmaxsize-jminsize, 0);
#endif
//printf("%d: DEBUG -- noffset %d cells %d\n",mype,noffset,ncells);
if (DEBUG) {
#ifdef _OPENMP
#pragma omp master
#endif
fprintf(fp,"%d: Sizes are imin %d imax %d jmin %d jmax %d\n",mype,iminsize,imaxsize,jminsize,jmaxsize);
}
static int imaxcalc, jmaxcalc;
#ifdef _OPENMP
#pragma omp for
#endif
for(uint ic=0; ic<ncells; ic++){
int cellnumber = ic+noffset;
int lev = level[ic];
int levmult = IPOW2(levmx-lev);
int ii = i[ic]*levmult-iminsize;
int jj = j[ic]*levmult-jminsize;
write_hash(cellnumber, jj*(imaxsize-iminsize)+ii, hash);
} // end for loop
if (TIMING_LEVEL >= 2) {
#ifdef _OPENMP
#pragma omp master
#endif
cpu_timers[MESH_TIMER_HASH_SETUP] += cpu_timer_stop(tstart_lev2);
cpu_timer_start(&tstart_lev2);
}
#ifdef _OPENMP
#pragma omp master
{
#endif
// Set neighbors to global cell numbers from hash
jmaxcalc = (jmax+1)*IPOW2(levmx);
imaxcalc = (imax+1)*IPOW2(levmx);
#ifdef _OPENMP
}
#pragma omp barrier
#endif
#ifdef _OPENMP
#pragma omp for
#endif
for (uint ic=0; ic<ncells; ic++){
int ii = i[ic];
int jj = j[ic];
int lev = level[ic];
int levmult = IPOW2(levmx-lev);
int iicur = ii*levmult-iminsize;
int iilft = max( (ii-1)*levmult, 0 )-iminsize;
int iirht = min( (ii+1)*levmult, imaxcalc-1)-iminsize;
int jjcur = jj*levmult-jminsize;
int jjbot = max( (jj-1)*levmult, 0 )-jminsize;
int jjtop = min( (jj+1)*levmult, jmaxcalc-1)-jminsize;
int nlftval = -1;
int nrhtval = -1;
int nbotval = -1;
int ntopval = -1;
// Taking care of boundary cells
// Force each boundary cell to point to itself on its boundary direction
if (iicur < 1*IPOW2(levmx) -iminsize) nlftval = ic+noffset;
if (jjcur < 1*IPOW2(levmx) -jminsize) nbotval = ic+noffset;
if (iicur > imax*IPOW2(levmx)-1-iminsize) nrhtval = ic+noffset;
if (jjcur > jmax*IPOW2(levmx)-1-jminsize) ntopval = ic+noffset;
// Boundary cells next to corner boundary need special checks
if (iicur == 1*IPOW2(levmx)-iminsize && (jjcur < 1*IPOW2(levmx)-jminsize || jjcur >= jmax*IPOW2(levmx)-jminsize ) ) nlftval = ic+noffset;
if (jjcur == 1*IPOW2(levmx)-jminsize && (iicur < 1*IPOW2(levmx)-iminsize || iicur >= imax*IPOW2(levmx)-iminsize ) ) nbotval = ic+noffset;
if (iirht == imax*IPOW2(levmx)-iminsize && (jjcur < 1*IPOW2(levmx)-jminsize || jjcur >= jmax*IPOW2(levmx)-jminsize ) ) nrhtval = ic+noffset;
if (jjtop == jmax*IPOW2(levmx)-jminsize && (iicur < 1*IPOW2(levmx)-iminsize || iicur >= imax*IPOW2(levmx)-iminsize ) ) ntopval = ic+noffset;
// need to check for finer neighbor first
// Right and top neighbor don't change for finer, so drop through to same size
// Left and bottom need to be half of same size index for finer test
if (lev != levmx) {
int iilftfiner = iicur-(iicur-iilft)/2;
int jjbotfiner = jjcur-(jjcur-jjbot)/2;
if (nlftval < 0) nlftval = read_hash(jjcur *(imaxsize-iminsize)+iilftfiner, hash);
if (nbotval < 0) nbotval = read_hash(jjbotfiner*(imaxsize-iminsize)+iicur, hash);
}
// same size neighbor
if (nlftval < 0) {
int nlfttry = read_hash(jjcur*(imaxsize-iminsize)+iilft, hash);
if (nlfttry >= 0 && nlfttry < (int)ncells && level[nlfttry] == lev) nlftval = nlfttry;
}
if (nrhtval < 0) nrhtval = read_hash(jjcur*(imaxsize-iminsize)+iirht, hash);
if (nbotval < 0) {
int nbottry = read_hash(jjbot*(imaxsize-iminsize)+iicur, hash);
if (nbottry >= 0 && nbottry < (int)ncells && level[nbottry] == lev) nbotval = nbottry;
}
if (ntopval < 0) ntopval = read_hash(jjtop*(imaxsize-iminsize)+iicur, hash);
// Now we need to take care of special case where bottom and left boundary need adjustment since
// expected cell doesn't exist on these boundaries if it is finer than current cell
if (lev != levmx) {
if (jjcur < 1*IPOW2(levmx)) {
if (nrhtval < 0) {
int jjtopfiner = (jjcur+jjtop)/2;
nrhtval = read_hash(jjtopfiner*(imaxsize-iminsize)+iirht, hash);
}
if (nlftval < 0) {
int iilftfiner = iicur-(iicur-iilft)/2;
int jjtopfiner = (jjcur+jjtop)/2;
nlftval = read_hash(jjtopfiner*(imaxsize-iminsize)+iilftfiner, hash);
}
}
if (iicur < 1*IPOW2(levmx)) {
if (ntopval < 0) {
int iirhtfiner = (iicur+iirht)/2;
ntopval = read_hash(jjtop*(imaxsize-iminsize)+iirhtfiner, hash);
}
if (nbotval < 0) {
int iirhtfiner = (iicur+iirht)/2;
int jjbotfiner = jjcur-(jjcur-jjbot)/2;
nbotval = read_hash(jjbotfiner*(imaxsize-iminsize)+iirhtfiner, hash);
}
}
}
// coarser neighbor
if (lev != 0){
if (nlftval < 0) {
iilft -= iicur-iilft;
int jjlft = (jj/2)*2*levmult-jminsize;
int nlfttry = read_hash(jjlft*(imaxsize-iminsize)+iilft, hash);
if (nlfttry >= 0 && nlfttry < (int)ncells && level[nlfttry] == lev-1) nlftval = nlfttry;
}
if (nrhtval < 0) {
int jjrht = (jj/2)*2*levmult-jminsize;
int nrhttry = read_hash(jjrht*(imaxsize-iminsize)+iirht, hash);
if (nrhttry >= 0 && nrhttry < (int)ncells && level[nrhttry] == lev-1) nrhtval = nrhttry;
}
if (nbotval < 0) {
jjbot -= jjcur-jjbot;
int iibot = (ii/2)*2*levmult-iminsize;
int nbottry = read_hash(jjbot*(imaxsize-iminsize)+iibot, hash);
if (nbottry >= 0 && nbottry < (int)ncells && level[nbottry] == lev-1) nbotval = nbottry;
}
if (ntopval < 0) {
int iitop = (ii/2)*2*levmult-iminsize;
int ntoptry = read_hash(jjtop*(imaxsize-iminsize)+iitop, hash);
if (ntoptry >= 0 && ntoptry < (int)ncells && level[ntoptry] == lev-1) ntopval = ntoptry;
}
}
nlft[ic] = nlftval;
nrht[ic] = nrhtval;
nbot[ic] = nbotval;
ntop[ic] = ntopval;
//fprintf(fp,"%d: neighbors[%d] = %d %d %d %d\n",mype,ic,nlft[ic],nrht[ic],nbot[ic],ntop[ic]);
}
if (DEBUG) {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
print_local();
int jmaxglobal = (jmax+1)*IPOW2(levmx);
int imaxglobal = (imax+1)*IPOW2(levmx);
fprintf(fp,"\n HASH 0 numbering\n");
for (int jj = jmaxglobal-1; jj>=0; jj--){
fprintf(fp,"%2d: %4d:",mype,jj);
if (jj >= jminsize && jj < jmaxsize) {
for (int ii = 0; ii<imaxglobal; ii++){
if (ii >= iminsize && ii < imaxsize) {
fprintf(fp,"%5d",read_hash((jj-jminsize)*(imaxsize-iminsize)+(ii-iminsize), hash));
} else {
fprintf(fp," ");
}
}
}
fprintf(fp,"\n");
}
fprintf(fp,"%2d: ",mype);
for (int ii = 0; ii<imaxglobal; ii++){
fprintf(fp,"%4d:",ii);
}
fprintf(fp,"\n");
fprintf(fp,"\n nlft numbering\n");
for (int jj = jmaxglobal-1; jj>=0; jj--){
fprintf(fp,"%2d: %4d:",mype,jj);
if (jj >= jminsize && jj < jmaxsize) {
for (int ii = 0; ii<imaxglobal; ii++){
if (ii >= iminsize && ii < imaxsize) {
int hashval = read_hash((jj-jminsize)*(imaxsize-iminsize)+(ii-iminsize), hash)-noffset;
if (hashval >= 0 && hashval < (int)ncells) {
fprintf(fp,"%5d",nlft[hashval]);
} else {
fprintf(fp," ");
}
} else {
fprintf(fp," ");
}
}
}
fprintf(fp,"\n");
}
fprintf(fp,"%2d: ",mype);
for (int ii = 0; ii<imaxglobal; ii++){
fprintf(fp,"%4d:",ii);
}
fprintf(fp,"\n");
fprintf(fp,"\n nrht numbering\n");
for (int jj = jmaxglobal-1; jj>=0; jj--){
fprintf(fp,"%2d: %4d:",mype,jj);
if (jj >= jminsize && jj < jmaxsize) {
for (int ii = 0; ii<imaxglobal; ii++){
if (ii >= iminsize && ii < imaxsize) {
int hashval = read_hash((jj-jminsize)*(imaxsize-iminsize)+(ii-iminsize), hash)-noffset;
if (hashval >= 0 && hashval < (int)ncells) {
fprintf(fp,"%5d",nrht[hashval]);
} else {
fprintf(fp," ");
}
} else {
fprintf(fp," ");
}
}
}
fprintf(fp,"\n");
}
fprintf(fp,"%2d: ",mype);
for (int ii = 0; ii<imaxglobal; ii++){
fprintf(fp,"%4d:",ii);
}
fprintf(fp,"\n");
fprintf(fp,"\n nbot numbering\n");
for (int jj = jmaxglobal-1; jj>=0; jj--){
fprintf(fp,"%2d: %4d:",mype,jj);
if (jj >= jminsize && jj < jmaxsize) {
for (int ii = 0; ii<imaxglobal; ii++){
if (ii >= iminsize && ii < imaxsize) {
int hashval = read_hash((jj-jminsize)*(imaxsize-iminsize)+(ii-iminsize), hash)-noffset;
if (hashval >= 0 && hashval < (int)ncells) {
fprintf(fp,"%5d",nbot[hashval]);
} else {
fprintf(fp," ");
}
} else {
fprintf(fp," ");
}
}
}
fprintf(fp,"\n");
}
fprintf(fp,"%2d: ",mype);
for (int ii = 0; ii<imaxglobal; ii++){
fprintf(fp,"%4d:",ii);
}
fprintf(fp,"\n");
fprintf(fp,"\n ntop numbering\n");
for (int jj = jmaxglobal-1; jj>=0; jj--){
fprintf(fp,"%2d: %4d:",mype,jj);
if (jj >= jminsize && jj < jmaxsize) {
for (int ii = 0; ii<imaxglobal; ii++){
if (ii >= iminsize && ii < imaxsize) {
int hashval = read_hash((jj-jminsize)*(imaxsize-iminsize)+(ii-iminsize), hash)-noffset;
if (hashval >= 0 && hashval < (int)ncells) {
fprintf(fp,"%5d",ntop[hashval]);
} else {
fprintf(fp," ");
}
} else {
fprintf(fp," ");
}
}
}
fprintf(fp,"\n");
}
fprintf(fp,"%2d: ",mype);
for (int ii = 0; ii<imaxglobal; ii++){
fprintf(fp,"%4d:",ii);
}
fprintf(fp,"\n");
#ifdef _OPENMP
}
#pragma omp barrier
#endif
}
if (TIMING_LEVEL >= 2) {
#ifdef _OPENMP
#pragma omp master
#endif
cpu_timers[MESH_TIMER_HASH_QUERY] += cpu_timer_stop(tstart_lev2);
cpu_timer_start(&tstart_lev2);
}
#ifdef HAVE_MPI
if (numpe > 1) {
static int num_comm_partners;
static vector<int> iminsize_global;
static vector<int> imaxsize_global;
static vector<int> jminsize_global;
static vector<int> jmaxsize_global;
static vector<int> comm_partner;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
iminsize_global.resize(numpe);
imaxsize_global.resize(numpe);
jminsize_global.resize(numpe);
jmaxsize_global.resize(numpe);
comm_partner.resize(numpe,-1);
MPI_Allgather(&iminsize, 1, MPI_INT, &iminsize_global[0], 1, MPI_INT, MPI_COMM_WORLD);
MPI_Allgather(&imaxsize, 1, MPI_INT, &imaxsize_global[0], 1, MPI_INT, MPI_COMM_WORLD);
MPI_Allgather(&jminsize, 1, MPI_INT, &jminsize_global[0], 1, MPI_INT, MPI_COMM_WORLD);
MPI_Allgather(&jmaxsize, 1, MPI_INT, &jmaxsize_global[0], 1, MPI_INT, MPI_COMM_WORLD);
num_comm_partners = 0;
for (int ip = 0; ip < numpe; ip++){
if (ip == mype) continue;
if (iminsize_global[ip] > imaxtile) continue;
if (imaxsize_global[ip] < imintile) continue;
if (jminsize_global[ip] > jmaxtile) continue;
if (jmaxsize_global[ip] < jmintile) continue;
comm_partner[num_comm_partners] = ip;
num_comm_partners++;
//if (DEBUG) fprintf(fp,"%d: overlap with processor %d bounding box is %d %d %d %d\n",mype,ip,iminsize_global[ip],imaxsize_global[ip],jminsize_global[ip],jmaxsize_global[ip]);
}
#ifdef _OPENMP
}
#pragma omp barrier
#endif
static vector<int> border_cell;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
border_cell.resize(ncells);
#ifdef BOUNDS_CHECK
for (uint ic=0; ic<ncells; ic++){
int nl = nlft[ic];
if (nl != -1){
nl -= noffset;
if (nl<0 || nl>= (int)ncells) printf("%d: Warning at line %d cell %d nlft %d\n",mype,__LINE__,ic,nl);
}
int nr = nrht[ic];
if (nr != -1){
nr -= noffset;
if (nr<0 || nr>= (int)ncells) printf("%d: Warning at line %d cell %d nrht %d\n",mype,__LINE__,ic,nr);
}
int nb = nbot[ic];
if (nb != -1){
nb -= noffset;
if (nb<0 || nb>= (int)ncells) printf("%d: Warning at line %d cell %d nbot %d\n",mype,__LINE__,ic,nb);
}
int nt = ntop[ic];
if (nt != -1){
nt -= noffset;
if (nt<0 || nt>= (int)ncells) printf("%d: Warning at line %d cell %d ntop %d\n",mype,__LINE__,ic,nt);
}
}
#endif
#ifdef _OPENMP
}
#pragma omp barrier
#endif
static vector<int> border_cell_out;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
border_cell_out.resize(ncells);
#ifdef _OPENMP
}
#pragma omp barrier
#endif
#ifdef _OPENMP
#pragma omp for
#endif
for (uint ic=0; ic<ncells; ic++){
int iborder_cell = 0;
// left neighbor is undefined -- or -- if left is at finer level check left top for undefined
if (nlft[ic] == -1 || (level[nlft[ic]-noffset] > level[ic] && ntop[nlft[ic]-noffset] == -1) ){
iborder_cell |= 0x0001;
}
if (nrht[ic] == -1 || (level[nrht[ic]-noffset] > level[ic] && ntop[nrht[ic]-noffset] == -1) ){
iborder_cell |= 0x0002;
}
if (nbot[ic] == -1 || (level[nbot[ic]-noffset] > level[ic] && nrht[nbot[ic]-noffset] == -1) ) {
iborder_cell |= 0x0004;
}
if (ntop[ic] == -1 || (level[ntop[ic]-noffset] > level[ic] && nrht[ntop[ic]-noffset] == -1) ) {
iborder_cell |= 0x0008;
}
border_cell[ic] = iborder_cell;
}
#ifdef _OPENMP
#pragma omp for
#endif
for (uint ic=0; ic<ncells; ic++){
int iborder_cell = border_cell[ic];
if (iborder_cell == 0) {
int nl = nlft[ic]-noffset;
if (nl >= 0 && nl < (int)ncells) {
if ((border_cell[nl] & 0x0001) == 0x0001) {
iborder_cell |= 0x0016;
} else if (level[nl] > level[ic]){
int ntl = ntop[nl]-noffset;
if (ntl >= 0 && ntl < (int)ncells && (border_cell[ntl] & 0x0001) == 0x0001) {
iborder_cell |= 0x0016;
}
}
}
int nr = nrht[ic]-noffset;
if (nr >= 0 && nr < (int)ncells) {
if ((border_cell[nrht[ic]-noffset] & 0x0002) == 0x0002) {
iborder_cell |= 0x0032;
} else if (level[nr] > level[ic]){
int ntr = ntop[nr]-noffset;
if (ntr >= 0 && ntr < (int)ncells && (border_cell[ntr] & 0x0002) == 0x0002) {
iborder_cell |= 0x0032;
}
}
}
int nb = nbot[ic]-noffset;
if (nb >= 0 && nb < (int)ncells) {
if ((border_cell[nb] & 0x0004) == 0x0004) {
iborder_cell |= 0x0064;
} else if (level[nb] > level[ic]){
int nrb = nrht[nb]-noffset;
if (nrb >= 0 && nrb < (int)ncells && (border_cell[nrb] & 0x0004) == 0x0004) {
iborder_cell |= 0x0064;
}
}
}
int nt = ntop[ic]-noffset;
if (nt >= 0 && nt < (int)ncells) {
if ((border_cell[nt] & 0x0008) == 0x0008) {
iborder_cell |= 0x0128;
} else if (level[nt] > level[ic]){
int nrt = nrht[nt]-noffset;
if (nrt >= 0 && nrt < (int)ncells && (border_cell[nrt] & 0x0008) == 0x0008) {
iborder_cell |= 0x0128;
}
}
}
}
border_cell_out[ic] = iborder_cell;
}
// indent offset
vector<int> border_cell_num;
static int nbsize_local;
static vector<int> border_cell_i;
static vector<int> border_cell_j;
static vector<int> border_cell_level;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
for (int ic=0; ic<(int)ncells; ic++){
if (border_cell_out[ic] > 0) border_cell_num.push_back(ic+noffset);
}
//printf("%d: border cell size is %d\n",mype,border_cell_num.size());
nbsize_local = border_cell_num.size();
border_cell_i.resize(nbsize_local);
border_cell_j.resize(nbsize_local);
border_cell_level.resize(nbsize_local);
for (int ic = 0; ic <nbsize_local; ic++){
int cell_num = border_cell_num[ic]-noffset;
border_cell_i[ic] = i[cell_num];
border_cell_j[ic] = j[cell_num];
border_cell_level[ic] = level[cell_num];
}
#ifdef _OPENMP
}
#pragma omp barrier
#endif
if (DEBUG) {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
fprintf(fp,"%d: Border cell size is %d\n",mype,nbsize_local);
for (int ib = 0; ib <nbsize_local; ib++){
fprintf(fp,"%d: Border cell %d is %d i %d j %d level %d\n",mype,ib,border_cell_num[ib],
border_cell_i[ib],border_cell_j[ib],border_cell_level[ib]);
}
#ifdef _OPENMP
}
#pragma omp barrier
#endif
}
if (TIMING_LEVEL >= 2) {
#ifdef _OPENMP
#pragma omp master
#endif
cpu_timers[MESH_TIMER_FIND_BOUNDARY] += cpu_timer_stop(tstart_lev2);
cpu_timer_start(&tstart_lev2);
}
// Allocate push database
static int **send_database;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
send_database = (int**)malloc(num_comm_partners*sizeof(int *));
for (int ip = 0; ip < num_comm_partners; ip++){
send_database[ip] = (int *)malloc(nbsize_local*sizeof(int));
}
#ifdef _OPENMP
}
#pragma omp barrier
#endif
// Compute the overlap between processor bounding boxes and set up push database
static vector<int> send_buffer_count;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
send_buffer_count.resize(num_comm_partners);
for (int ip = 0; ip < num_comm_partners; ip++){
int icount = 0;
for (int ib = 0; ib <nbsize_local; ib++){
int lev = border_cell_level[ib];
int levmult = IPOW2(levmx-lev);
if (border_cell_i[ib]*levmult >= iminsize_global[comm_partner[ip]] &&
border_cell_i[ib]*levmult <= imaxsize_global[comm_partner[ip]] &&
border_cell_j[ib]*levmult >= jminsize_global[comm_partner[ip]] &&
border_cell_j[ib]*levmult <= jmaxsize_global[comm_partner[ip]] ) {
// border_cell_i[ib],border_cell_j[ib],border_cell_level[ib]);
send_database[ip][icount] = ib;
icount++;
}
}
send_buffer_count[ip]=icount;
}
#ifdef _OPENMP
}
#pragma omp barrier
#endif
// Initialize L7_Push_Setup with num_comm_partners, comm_partner, send_database and
// send_buffer_count. L7_Push_Setup will copy data and determine recv_buffer_counts.
// It will return receive_count_total for use in allocations
static int receive_count_total;
int i_push_handle = 0;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
i_push_handle = 0;
L7_Push_Setup(num_comm_partners, &comm_partner[0], &send_buffer_count[0],
send_database, &receive_count_total, &i_push_handle);
#ifdef _OPENMP
}
#pragma omp barrier
#endif
if (DEBUG) {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
fprintf(fp,"DEBUG num_comm_partners %d\n",num_comm_partners);
for (int ip = 0; ip < num_comm_partners; ip++){
fprintf(fp,"DEBUG comm partner is %d data count is %d\n",comm_partner[ip],send_buffer_count[ip]);
for (int ic = 0; ic < send_buffer_count[ip]; ic++){
int ib = send_database[ip][ic];
fprintf(fp,"DEBUG \t index %d cell number %d i %d j %d level %d\n",ib,border_cell_num[ib],
border_cell_i[ib],border_cell_j[ib],border_cell_level[ib]);
}
}
#ifdef _OPENMP
}
#endif
}
// Can now free the send database. Other arrays are vectors and will automatically
// deallocate
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
for (int ip = 0; ip < num_comm_partners; ip++){
free(send_database[ip]);
}
free(send_database);
#ifdef _OPENMP
}
#endif
if (TIMING_LEVEL >= 2) {
#ifdef _OPENMP
#pragma omp master
#endif
cpu_timers[MESH_TIMER_PUSH_SETUP] += cpu_timer_stop(tstart_lev2);
cpu_timer_start(&tstart_lev2);
}
// Push the data needed to the adjacent processors
static int *border_cell_num_local;
static int *border_cell_i_local;
static int *border_cell_j_local;
static int *border_cell_level_local;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
border_cell_num_local = (int *)malloc(receive_count_total*sizeof(int));
border_cell_i_local = (int *)malloc(receive_count_total*sizeof(int));
border_cell_j_local = (int *)malloc(receive_count_total*sizeof(int));
border_cell_level_local = (int *)malloc(receive_count_total*sizeof(int));
L7_Push_Update(&border_cell_num[0], border_cell_num_local, i_push_handle);
L7_Push_Update(&border_cell_i[0], border_cell_i_local, i_push_handle);
L7_Push_Update(&border_cell_j[0], border_cell_j_local, i_push_handle);
L7_Push_Update(&border_cell_level[0], border_cell_level_local, i_push_handle);
L7_Push_Free(&i_push_handle);
#ifdef _OPENMP
}
#pragma omp barrier
#endif
nbsize_local = receive_count_total;
if (DEBUG) {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
for (int ic = 0; ic < nbsize_local; ic++) {
fprintf(fp,"%d: Local Border cell %d is %d i %d j %d level %d\n",mype,ic,border_cell_num_local[ic],
border_cell_i_local[ic],border_cell_j_local[ic],border_cell_level_local[ic]);
}
#ifdef _OPENMP
}
#endif
}
if (TIMING_LEVEL >= 2) {
#ifdef _OPENMP
#pragma omp master
#endif
cpu_timers[MESH_TIMER_PUSH_BOUNDARY] += cpu_timer_stop(tstart_lev2);
cpu_timer_start(&tstart_lev2);
}
if (TIMING_LEVEL >= 2) {
#ifdef _OPENMP
#pragma omp master
#endif
cpu_timers[MESH_TIMER_LOCAL_LIST] += cpu_timer_stop(tstart_lev2);
cpu_timer_start(&tstart_lev2);
}
if (DEBUG) {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
int jmaxglobal = (jmax+1)*IPOW2(levmx);
int imaxglobal = (imax+1)*IPOW2(levmx);
fprintf(fp,"\n HASH numbering before layer 1\n");
for (int jj = jmaxglobal-1; jj>=0; jj--){
fprintf(fp,"%2d: %4d:",mype,jj);
if (jj >= jminsize && jj < jmaxsize) {
for (int ii = 0; ii<imaxglobal; ii++){
if (ii >= iminsize && ii < imaxsize) {
fprintf(fp,"%5d",read_hash((jj-jminsize)*(imaxsize-iminsize)+(ii-iminsize), hash));
} else {
fprintf(fp," ");
}
}
}
fprintf(fp,"\n");
}
fprintf(fp,"%2d: ",mype);
for (int ii = 0; ii<imaxglobal; ii++){
fprintf(fp,"%4d:",ii);
}
fprintf(fp,"\n");
#ifdef _OPENMP
}
#endif
}
vector<int> border_cell_needed_local;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
border_cell_needed_local.resize(nbsize_local, 0);
#ifdef _OPENMP
}
#pragma omp barrier
#endif
// Layer 1
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
for (int ic =0; ic<nbsize_local; ic++){
int jj = border_cell_j_local[ic];
int ii = border_cell_i_local[ic];
int lev = border_cell_level_local[ic];
int levmult = IPOW2(levmx-lev);
int iicur = ii*levmult-iminsize;
int iilft = max( (ii-1)*levmult, 0 )-iminsize;
int iirht = min( (ii+1)*levmult, imaxcalc-1)-iminsize;
int jjcur = jj*levmult-jminsize;
int jjbot = max( (jj-1)*levmult, 0 )-jminsize;
int jjtop = min( (jj+1)*levmult, jmaxcalc-1)-jminsize;
//fprintf(fp,"DEBUG layer ic %d num %d i %d j %d lev %d\n",ic,border_cell_num_local[ic],ii,jj,lev);
int iborder = 0;
// Test for cell to left
if (iicur-(iicur-iilft)/2 >= 0 && iicur-(iicur-iilft)/2 < imaxsize-iminsize && jjcur >= 0 && (jjcur+jjtop)/2 < jmaxsize-jminsize){
int nlftval = -1;
// Check for finer cell left and bottom side
if (lev != levmx){ // finer neighbor
int iilftfiner = iicur-(iicur-iilft)/2;
nlftval = read_hash(jjcur*(imaxsize-iminsize)+iilftfiner, hash);
// Also check for finer cell left and top side
if (nlftval < 0) {
int jjtopfiner = (jjcur+jjtop)/2;
nlftval = read_hash(jjtopfiner*(imaxsize-iminsize)+iilftfiner, hash);
}
}
if (nlftval < 0 && iilft >= 0) { // same size
int nlfttry = read_hash(jjcur*(imaxsize-iminsize)+iilft, hash);
// we have to test for same level or it could be a finer cell one cell away that it is matching
if (nlfttry-noffset >= 0 && nlfttry-noffset < (int)ncells && level[nlfttry-noffset] == lev) {
nlftval = nlfttry;
}
}
if (lev != 0 && nlftval < 0 && iilft-(iicur-iilft) >= 0){ // coarser neighbor
iilft -= iicur-iilft;
int jjlft = (jj/2)*2*levmult-jminsize;
int nlfttry = read_hash(jjlft*(imaxsize-iminsize)+iilft, hash);
// we have to test for coarser level or it could be a same size cell one or two cells away that it is matching
if (nlfttry-noffset >= 0 && nlfttry-noffset < (int)ncells && level[nlfttry-noffset] == lev-1) {
nlftval = nlfttry;
}
}
if (nlftval >= 0) iborder |= 0x0001;
}
// Test for cell to right
if (iirht < imaxsize-iminsize && iirht >= 0 && jjcur >= 0 && jjtop < jmaxsize-jminsize) {
int nrhtval = -1;
// right neighbor -- finer, same size and coarser
nrhtval = read_hash(jjcur*(imaxsize-iminsize)+iirht, hash);
// right neighbor -- finer right top test
if (nrhtval < 0 && lev != levmx){
int jjtopfiner = (jjcur+jjtop)/2;
nrhtval = read_hash(jjtopfiner*(imaxsize-iminsize)+iirht, hash);
}
if (nrhtval < 0 && lev != 0) { // test for coarser, but not directly above
int jjrhtcoarser = (jj/2)*2*levmult-jminsize;
if (jjrhtcoarser != jjcur) {
int nrhttry = read_hash(jjrhtcoarser*(imaxsize-iminsize)+iirht, hash);
if (nrhttry-noffset >= 0 && nrhttry-noffset < (int)ncells && level[nrhttry-noffset] == lev-1) {
nrhtval = nrhttry;
}
}
}
if (nrhtval > 0) iborder |= 0x0002;
}
// Test for cell to bottom
if (iicur >= 0 && (iicur+iirht)/2 < imaxsize-iminsize && jjcur-(jjcur-jjbot)/2 >= 0 && jjcur-(jjcur-jjbot)/2 < jmaxsize-jminsize){
int nbotval = -1;
// Check for finer cell below and left side
if (lev != levmx){ // finer neighbor
int jjbotfiner = jjcur-(jjcur-jjbot)/2;
nbotval = read_hash(jjbotfiner*(imaxsize-iminsize)+iicur, hash);
// Also check for finer cell below and right side
if (nbotval < 0) {
int iirhtfiner = (iicur+iirht)/2;
nbotval = read_hash(jjbotfiner*(imaxsize-iminsize)+iirhtfiner, hash);
}
}
if (nbotval < 0 && jjbot >= 0) { // same size
int nbottry = read_hash(jjbot*(imaxsize-iminsize)+iicur, hash);
// we have to test for same level or it could be a finer cell one cell away that it is matching
if (nbottry-noffset >= 0 && nbottry-noffset < (int)ncells && level[nbottry-noffset] == lev) {
nbotval = nbottry;
}
}
if (lev != 0 && nbotval < 0 && jjbot-(jjcur-jjbot) >= 0){ // coarser neighbor
jjbot -= jjcur-jjbot;
int iibot = (ii/2)*2*levmult-iminsize;
int nbottry = read_hash(jjbot*(imaxsize-iminsize)+iibot, hash);
// we have to test for coarser level or it could be a same size cell one or two cells away that it is matching
if (nbottry-noffset >= 0 && nbottry-noffset < (int)ncells && level[nbottry-noffset] == lev-1) {
nbotval = nbottry;
}
}
if (nbotval >= 0) iborder |= 0x0004;
}
// Test for cell to top
if (iirht < imaxsize-iminsize && iicur >= 0 && jjtop >= 0 && jjtop < jmaxsize-jminsize) {
int ntopval = -1;
// top neighbor -- finer, same size and coarser
ntopval = read_hash(jjtop*(imaxsize-iminsize)+iicur, hash);
// top neighbor -- finer top right test
if (ntopval < 0 && lev != levmx){
int iirhtfiner = (iicur+iirht)/2;
ntopval = read_hash(jjtop*(imaxsize-iminsize)+iirhtfiner, hash);
}
if (ntopval < 0 && lev != 0) { // test for coarser, but not directly above
int iitopcoarser = (ii/2)*2*levmult-iminsize;
if (iitopcoarser != iicur) {
int ntoptry = read_hash(jjtop*(imaxsize-iminsize)+iitopcoarser, hash);
if (ntoptry-noffset >= 0 && ntoptry-noffset < (int)ncells && level[ntoptry-noffset] == lev-1) {
ntopval = ntoptry;
}
}
}
if (ntopval > 0) iborder |= 0x0008;
}
if (iborder) border_cell_needed_local[ic] = iborder;
}
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
if (DEBUG) {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
for(int ic=0; ic<nbsize_local; ic++){
if (border_cell_needed_local[ic] == 0) continue;
fprintf(fp,"%d: First set of needed cells ic %3d cell %3d type %3d\n",mype,ic,border_cell_num_local[ic],border_cell_needed_local[ic]);
}
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
}
// Walk through cell array and set hash to border local index plus ncells+noffset for next pass
//fprintf(fp,"%d: DEBUG new hash jminsize %d jmaxsize %d iminsize %d imaxsize %d\n",mype,jminsize,jmaxsize,iminsize,imaxsize);
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
for(int ic=0; ic<nbsize_local; ic++){
if (border_cell_needed_local[ic] == 0) continue;
//fprintf(fp,"%d: index %d cell %d i %d j %d\n",mype,ic,border_cell_num_local[ic],border_cell_i_local[ic],border_cell_j_local[ic]);
int lev = border_cell_level_local[ic];
int levmult = IPOW2(levmx-lev);
int ii = border_cell_i_local[ic]*levmult-iminsize;
int jj = border_cell_j_local[ic]*levmult-jminsize;
write_hash(ncells+noffset+ic, jj*(imaxsize-iminsize)+ii, hash);
}
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
if (TIMING_LEVEL >= 2) {
#ifdef _OPENMP
#pragma omp master
#endif
cpu_timers[MESH_TIMER_LAYER1] += cpu_timer_stop(tstart_lev2);
cpu_timer_start(&tstart_lev2);
}
if (DEBUG) {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
print_local();
int jmaxglobal = (jmax+1)*IPOW2(levmx);
int imaxglobal = (imax+1)*IPOW2(levmx);
fprintf(fp,"\n HASH numbering for 1 layer\n");
for (int jj = jmaxglobal-1; jj>=0; jj--){
fprintf(fp,"%2d: %4d:",mype,jj);
if (jj >= jminsize && jj < jmaxsize) {
for (int ii = 0; ii<imaxglobal; ii++){
if (ii >= iminsize && ii < imaxsize) {
fprintf(fp,"%5d",read_hash((jj-jminsize)*(imaxsize-iminsize)+(ii-iminsize), hash) );
} else {
fprintf(fp," ");
}
}
}
fprintf(fp,"\n");
}
fprintf(fp,"%2d: ",mype);
for (int ii = 0; ii<imaxglobal; ii++){
fprintf(fp,"%4d:",ii);
}
fprintf(fp,"\n");
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
}
// Layer 2
#ifdef _OPENMP
#pragma omp master
{
#endif
for (int ic =0; ic<nbsize_local; ic++){
if (border_cell_needed_local[ic] > 0) continue;
int jj = border_cell_j_local[ic];
int ii = border_cell_i_local[ic];
int lev = border_cell_level_local[ic];
int levmult = IPOW2(levmx-lev);
int iicur = ii*levmult-iminsize;
int iilft = max( (ii-1)*levmult, 0 )-iminsize;
int iirht = min( (ii+1)*levmult, imaxcalc-1)-iminsize;
int jjcur = jj*levmult-jminsize;
int jjbot = max( (jj-1)*levmult, 0 )-jminsize;
int jjtop = min( (jj+1)*levmult, jmaxcalc-1)-jminsize;
//fprintf(fp," DEBUG layer2 ic %d num %d i %d j %d lev %d\n",ic,border_cell_num_local[ic],ii,jj,lev);
int iborder = 0;
// Test for cell to left
if (iicur-(iicur-iilft)/2 >= 0 && iicur-(iicur-iilft)/2 < imaxsize-iminsize && jjcur >= 0 && (jjcur+jjtop)/2 < jmaxsize-jminsize){
// Check for finer cell left and bottom side
if (lev != levmx){ // finer neighbor
int iilftfiner = iicur-(iicur-iilft)/2;
int nl = read_hash(jjcur*(imaxsize-iminsize)+iilftfiner, hash);
if (nl >= (int)(ncells+noffset) && (border_cell_needed_local[nl-ncells-noffset] & 0x0001) == 0x0001) {
iborder = 0x0001;
} else {
// Also check for finer cell left and top side
int jjtopfiner = (jjcur+jjtop)/2;
int nlt = read_hash(jjtopfiner*(imaxsize-iminsize)+iilftfiner, hash);
if ( nlt >= (int)(ncells+noffset) && (border_cell_needed_local[nlt-ncells-noffset] & 0x0001) == 0x0001) {
iborder = 0x0001;
}
}
}
if ( (iborder & 0x0001) == 0 && iilft >= 0) { //same size
int nl = read_hash(jjcur*(imaxsize-iminsize)+iilft, hash);
int levcheck = -1;
if (nl-noffset >= 0 && nl-noffset < (int)ncells) {
levcheck = level[nl-noffset];
} else if (nl >= 0 && (int)(nl-ncells-noffset) >= 0 && (int)(nl-ncells-noffset) < nbsize_local) {
levcheck = border_cell_level_local[nl-ncells-noffset];
}
if (nl >= (int)(ncells+noffset) && levcheck == lev && (border_cell_needed_local[nl-ncells-noffset] & 0x0001) == 0x0001) {
iborder = 0x0001;
} else if (lev != 0 && iilft-(iicur-iilft) >= 0){ // coarser neighbor
iilft -= iicur-iilft;
int jjlft = (jj/2)*2*levmult-jminsize;
nl = read_hash(jjlft*(imaxsize-iminsize)+iilft, hash);
levcheck = -1;
if (nl-noffset >= 0 && nl-noffset < (int)ncells) {
levcheck = level[nl-noffset];
} else if (nl >= 0 && (int)(nl-ncells-noffset) >= 0 && (int)(nl-ncells-noffset) < nbsize_local) {
levcheck = border_cell_level_local[nl-ncells-noffset];
}
// we have to test for coarser level or it could be a same size cell one or two cells away that it is matching
if (nl >= (int)(ncells+noffset) && levcheck == lev-1 && (border_cell_needed_local[nl-ncells-noffset] & 0x0001) == 0x0001) {
iborder = 0x0001;
}
}
}
}
// Test for cell to right
if (iirht < imaxsize-iminsize && iirht >= 0 && jjcur >= 0 && jjtop < jmaxsize-jminsize) {
// right neighbor -- finer, same size and coarser
int nr = read_hash(jjcur*(imaxsize-iminsize)+iirht, hash);
if (nr >= (int)(ncells+noffset) && (border_cell_needed_local[nr-ncells-noffset] & 0x0002) == 0x0002) {
iborder = 0x0002;
} else if (lev != levmx){
// right neighbor -- finer right top test
int jjtopfiner = (jjcur+jjtop)/2;
int nrt = read_hash(jjtopfiner*(imaxsize-iminsize)+iirht, hash);
if (nrt >= (int)(ncells+noffset) && (border_cell_needed_local[nrt-ncells-noffset] & 0x0002) == 0x0002) {
iborder = 0x0002;
}
}
if ( (iborder & 0x0002) == 0 && lev != 0) { // test for coarser, but not directly right
int jjrhtcoarser = (jj/2)*2*levmult-jminsize;
if (jjrhtcoarser != jjcur) {
int nr = read_hash(jjrhtcoarser*(imaxsize-iminsize)+iirht, hash);
int levcheck = -1;
if (nr-noffset >= 0 && nr-noffset < (int)ncells) {
levcheck = level[nr-noffset];
} else if (nr >= 0 && (int)(nr-ncells-noffset) >= 0 && (int)(nr-ncells-noffset) < nbsize_local) {
levcheck = border_cell_level_local[nr-ncells-noffset];
}
if (nr >= (int)(ncells+noffset) && levcheck == lev-1 && (border_cell_needed_local[nr-ncells-noffset] & 0x0002) == 0x0002) {
iborder = 0x0002;
}
}
}
}
// Test for cell to bottom
if (iicur >= 0 && (iicur+iirht)/2 < imaxsize-iminsize && jjcur-(jjcur-jjbot)/2 >= 0 && jjcur-(jjcur-jjbot)/2 < jmaxsize-jminsize){
// Check for finer cell below and left side
if (lev != levmx){ // finer neighbor
int jjbotfiner = jjcur-(jjcur-jjbot)/2;
int nb = read_hash(jjbotfiner*(imaxsize-iminsize)+iicur, hash);
if (nb >= (int)(ncells+noffset) && (border_cell_needed_local[nb-ncells-noffset] & 0x0004) == 0x0004) {
iborder = 0x0004;
} else {
// Also check for finer cell below and right side
int iirhtfiner = (iicur+iirht)/2;
int nbr = read_hash(jjbotfiner*(imaxsize-iminsize)+iirhtfiner, hash);
if (nbr >= (int)(ncells+noffset) && (border_cell_needed_local[nbr-ncells-noffset] & 0x0004) == 0x0004) {
iborder = 0x0004;
}
}
}
if ( (iborder & 0x0004) == 0 && jjbot >= 0) { //same size
int nb = read_hash(jjbot*(imaxsize-iminsize)+iicur, hash);
int levcheck = -1;
if (nb-noffset >= 0 && nb-noffset < (int)ncells) {
levcheck = level[nb-noffset];
} else if (nb >= 0 && (int)(nb-ncells-noffset) >= 0 && (int)(nb-ncells-noffset) < nbsize_local) {
levcheck = border_cell_level_local[nb-ncells-noffset];
}
if (nb >= (int)(ncells+noffset) && levcheck == lev && (border_cell_needed_local[nb-ncells-noffset] & 0x0004) == 0x0004) {
iborder = 0x0004;
} else if (lev != 0 && jjbot-(jjcur-jjbot) >= 0){ // coarser neighbor
jjbot -= jjcur-jjbot;
int iibot = (ii/2)*2*levmult-iminsize;
nb = read_hash(jjbot*(imaxsize-iminsize)+iibot, hash);
levcheck = -1;
if (nb-noffset >= 0 && nb-noffset < (int)ncells) {
levcheck = level[nb-noffset];
} else if (nb >= 0 && (int)(nb-ncells-noffset) >= 0 && (int)(nb-ncells-noffset) < nbsize_local) {
levcheck = border_cell_level_local[nb-ncells-noffset];
}
// we have to test for coarser level or it could be a same size cell one or two cells away that it is matching
if (nb >= (int)(ncells+noffset) && levcheck == lev-1 && (border_cell_needed_local[nb-ncells-noffset] & 0x0004) == 0x0004) {
iborder = 0x0004;
}
}
}
}
// Test for cell to top
if (iirht < imaxsize-iminsize && iicur >= 0 && jjtop >= 0 && jjtop < jmaxsize-jminsize) {
// top neighbor -- finer, same size and coarser
int nt = read_hash(jjtop*(imaxsize-iminsize)+iicur, hash);
if (nt >= (int)(ncells+noffset) && (border_cell_needed_local[nt-ncells-noffset] & 0x0008) == 0x0008) {
iborder = 0x0008;
} else if (lev != levmx){
int iirhtfiner = (iicur+iirht)/2;
int ntr = read_hash(jjtop*(imaxsize-iminsize)+iirhtfiner, hash);
if ( ntr >= (int)(ncells+noffset) && (border_cell_needed_local[ntr-ncells-noffset] & 0x0008) == 0x0008) {
iborder = 0x0008;
}
}
if ( (iborder & 0x0008) == 0 && lev != 0) { // test for coarser, but not directly above
int iitopcoarser = (ii/2)*2*levmult-iminsize;
if (iitopcoarser != iicur) {
int nb = read_hash(jjtop*(imaxsize-iminsize)+iitopcoarser, hash);
int levcheck = -1;
if (nb-noffset >= 0 && nb-noffset < (int)ncells) {
levcheck = level[nb-noffset];
} else if (nb >= 0 && (int)(nb-ncells-noffset) >= 0 && (int)(nb-ncells-noffset) < nbsize_local) {
levcheck = border_cell_level_local[nb-ncells-noffset];
}
if (nb-noffset >= (int)(ncells-noffset) && levcheck == lev-1 && (border_cell_needed_local[nb-ncells-noffset] & 0x0008) == 0x0008) {
iborder = 0x0008;
}
}
}
}
if (iborder) border_cell_needed_local[ic] = iborder |= 0x0016;
}
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
vector<int> indices_needed;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
if (DEBUG) {
for(int ic=0; ic<nbsize_local; ic++){
if (border_cell_needed_local[ic] < 0x0016) fprintf(fp,"%d: First set of needed cells ic %3d cell %3d type %3d\n",mype,ic,border_cell_num_local[ic],border_cell_needed_local[ic]);
if (border_cell_needed_local[ic] >= 0x0016) fprintf(fp,"%d: Second set of needed cells ic %3d cell %3d type %3d\n",mype,ic,border_cell_num_local[ic],border_cell_needed_local[ic]);
}
}
int inew = 0;
for(int ic=0; ic<nbsize_local; ic++){
if (border_cell_needed_local[ic] <= 0) continue;
indices_needed.push_back(border_cell_num_local[ic]);
border_cell_num_local[inew] = border_cell_num_local[ic];
border_cell_i_local[inew] = border_cell_i_local[ic];
border_cell_j_local[inew] = border_cell_j_local[ic];
border_cell_level_local[inew] = border_cell_level_local[ic];
// border_cell_num_local is not used after -- could be commented out?
// border_cell_needed_local[inew] = 1;
inew++;
}
nbsize_local = inew;
free(border_cell_num_local);
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
// Walk through cell array and set hash to global cell values
//fprintf(fp,"%d: DEBUG new hash jminsize %d jmaxsize %d iminsize %d imaxsize %d\n",mype,jminsize,jmaxsize,iminsize,imaxsize);
#ifdef _OPENMP
#pragma omp for
#endif
for(int ic=0; ic<nbsize_local; ic++){
int lev = border_cell_level_local[ic];
int levmult = IPOW2(levmx-lev);
int ii = border_cell_i_local[ic]*levmult-iminsize;
int jj = border_cell_j_local[ic]*levmult-jminsize;
write_hash(-(ncells+ic), jj*(imaxsize-iminsize)+ii, hash);
}
if (TIMING_LEVEL >= 2) {
#ifdef _OPENMP
#pragma omp master
#endif
cpu_timers[MESH_TIMER_LAYER2] += cpu_timer_stop(tstart_lev2);
cpu_timer_start(&tstart_lev2);
}
if (DEBUG) {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
print_local();
int jmaxglobal = (jmax+1)*IPOW2(levmx);
int imaxglobal = (imax+1)*IPOW2(levmx);
fprintf(fp,"\n HASH numbering for 2 layer\n");
for (int jj = jmaxglobal-1; jj>=0; jj--){
fprintf(fp,"%2d: %4d:",mype,jj);
if (jj >= jminsize && jj < jmaxsize) {
for (int ii = 0; ii<imaxglobal; ii++){
if (ii >= iminsize && ii < imaxsize) {
fprintf(fp,"%5d",read_hash((jj-jminsize)*(imaxsize-iminsize)+(ii-iminsize), hash) );
} else {
fprintf(fp," ");
}
}
}
fprintf(fp,"\n");
}
fprintf(fp,"%2d: ",mype);
for (int ii = 0; ii<imaxglobal; ii++){
fprintf(fp,"%4d:",ii);
}
fprintf(fp,"\n");
#ifdef _OPENMP
} // end master region
#endif
}
if (TIMING_LEVEL >= 2) {
#ifdef _OPENMP
#pragma omp master
#endif
cpu_timers[MESH_TIMER_LAYER_LIST] += cpu_timer_stop(tstart_lev2);
cpu_timer_start(&tstart_lev2);
}
int nghost = nbsize_local;
ncells_ghost = ncells + nghost;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
celltype = (int *)mesh_memory.memory_realloc(ncells_ghost, celltype);
i = (int *)mesh_memory.memory_realloc(ncells_ghost, i);
j = (int *)mesh_memory.memory_realloc(ncells_ghost, j);
level = (int *)mesh_memory.memory_realloc(ncells_ghost, level);
nlft = (int *)mesh_memory.memory_realloc(ncells_ghost, nlft);
nrht = (int *)mesh_memory.memory_realloc(ncells_ghost, nrht);
nbot = (int *)mesh_memory.memory_realloc(ncells_ghost, nbot);
ntop = (int *)mesh_memory.memory_realloc(ncells_ghost, ntop);
memory_reset_ptrs();
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
#ifdef _OPENMP
#pragma omp for
#endif
for (int ic = ncells; ic < (int)ncells_ghost; ic++){
nlft[ic] = -1;
nrht[ic] = -1;
nbot[ic] = -1;
ntop[ic] = -1;
}
if (TIMING_LEVEL >= 2) {
#ifdef _OPENMP
#pragma omp master
#endif
cpu_timers[MESH_TIMER_COPY_MESH_DATA] += cpu_timer_stop(tstart_lev2);
cpu_timer_start(&tstart_lev2);
}
#ifdef _OPENMP
#pragma omp for
#endif
for(int ic=0; ic<nbsize_local; ic++){
int ii = border_cell_i_local[ic];
int jj = border_cell_j_local[ic];
int lev = border_cell_level_local[ic];
if (ii < lev_ibegin[lev]) celltype[ncells+ic] = LEFT_BOUNDARY;
if (ii > lev_iend[lev]) celltype[ncells+ic] = RIGHT_BOUNDARY;
if (jj < lev_jbegin[lev]) celltype[ncells+ic] = BOTTOM_BOUNDARY;
if (jj > lev_jend[lev]) celltype[ncells+ic] = TOP_BOUNDARY;
i[ncells+ic] = ii;
j[ncells+ic] = jj;
level[ncells+ic] = lev;
}
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
free(border_cell_i_local);
free(border_cell_j_local);
free(border_cell_level_local);
#ifdef _OPENMP
} // end master region
#endif
if (TIMING_LEVEL >= 2) {
#ifdef _OPENMP
#pragma omp master
#endif
cpu_timers[MESH_TIMER_FILL_MESH_GHOST] += cpu_timer_stop(tstart_lev2);
cpu_timer_start(&tstart_lev2);
}
if (DEBUG) {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
fprintf(fp,"After copying i,j, level to ghost cells\n");
print_local();
#ifdef _OPENMP
} // end master region
#endif
}
#ifdef _OPENMP
#pragma omp for
#endif
for (uint ic=0; ic<ncells_ghost; ic++){
int ii = i[ic];
int jj = j[ic];
int lev = level[ic];
int levmult = IPOW2(levmx-lev);
int iicur = ii*levmult-iminsize;
int iilft = max( (ii-1)*levmult, 0 )-iminsize;
int iirht = min( (ii+1)*levmult, imaxcalc-1)-iminsize;
int jjcur = jj*levmult-jminsize;
int jjbot = max( (jj-1)*levmult, 0 )-jminsize;
int jjtop = min( (jj+1)*levmult, jmaxcalc-1)-jminsize;
//fprintf(fp,"DEBUG neigh ic %d nlft %d ii %d levmult %d iminsize %d icheck %d\n",ic,nlft[ic],ii,levmult,iminsize,(max( ii *levmult-1, 0))-iminsize);
int nlftval = nlft[ic];
int nrhtval = nrht[ic];
int nbotval = nbot[ic];
int ntopval = ntop[ic];
if (nlftval == -1){
// Taking care of boundary cells
// Force each boundary cell to point to itself on its boundary direction
if (iicur < 1*IPOW2(levmx) -iminsize) nlftval = read_hash(jjcur*(imaxsize-iminsize)+iicur, hash);
// Boundary cells next to corner boundary need special checks
if (iicur == 1*IPOW2(levmx)-iminsize && (jjcur < 1*IPOW2(levmx)-jminsize || jjcur >= jmax*IPOW2(levmx)-jminsize ) ) nlftval = read_hash(jjcur*(imaxsize-iminsize)+iicur, hash);
// need to check for finer neighbor first
// Right and top neighbor don't change for finer, so drop through to same size
// Left and bottom need to be half of same size index for finer test
if (lev != levmx) {
int iilftfiner = iicur-(iicur-iilft)/2;
if (nlftval == -1 && iilftfiner >= 0) nlftval = read_hash(jjcur*(imaxsize-iminsize)+iilftfiner, hash);
}
// same size neighbor
if (nlftval == -1 && iilft >= 0) nlftval = read_hash(jjcur*(imaxsize-iminsize)+iilft, hash);
// Now we need to take care of special case where bottom and left boundary need adjustment since
// expected cell doesn't exist on these boundaries if it is finer than current cell
if (jjcur < 1*IPOW2(levmx) && lev != levmx) {
if (nlftval == -1) {
int iilftfiner = iicur-(iicur-iilft)/2;
int jjtopfiner = (jjcur+jjtop)/2;
if (jjtopfiner < jmaxsize-jminsize && iilftfiner >= 0) nlftval = read_hash(jjtopfiner*(imaxsize-iminsize)+iilftfiner, hash);
}
}
// coarser neighbor
if (lev != 0){
if (nlftval == -1) {
int iilftcoarser = iilft - (iicur-iilft);
int jjlft = (jj/2)*2*levmult-jminsize;
if (iilftcoarser >=0) nlftval = read_hash(jjlft*(imaxsize-iminsize)+iilftcoarser, hash);
}
}
if (nlftval != -1) nlft[ic] = nlftval;
}
if (nrhtval == -1) {
// Taking care of boundary cells
// Force each boundary cell to point to itself on its boundary direction
if (iicur > imax*IPOW2(levmx)-1-iminsize) nrhtval = read_hash(jjcur*(imaxsize-iminsize)+iicur, hash);
// Boundary cells next to corner boundary need special checks
if (iirht == imax*IPOW2(levmx)-iminsize && (jjcur < 1*IPOW2(levmx)-jminsize || jjcur >= jmax*IPOW2(levmx)-jminsize ) ) nrhtval = read_hash(jjcur*(imaxsize-iminsize)+iicur, hash);
// same size neighbor
if (nrhtval == -1 && iirht < imaxsize-iminsize) nrhtval = read_hash(jjcur*(imaxsize-iminsize)+iirht, hash);
// Now we need to take care of special case where bottom and left boundary need adjustment since
// expected cell doesn't exist on these boundaries if it is finer than current cell
if (jjcur < 1*IPOW2(levmx) && lev != levmx) {
if (nrhtval == -1) {
int jjtopfiner = (jjcur+jjtop)/2;
if (jjtopfiner < jmaxsize-jminsize && iirht < imaxsize-iminsize) nrhtval = read_hash(jjtopfiner*(imaxsize-iminsize)+iirht, hash);
}
}
// coarser neighbor
if (lev != 0){
if (nrhtval == -1) {
int jjrht = (jj/2)*2*levmult-jminsize;
if (iirht < imaxsize-iminsize) nrhtval = read_hash(jjrht*(imaxsize-iminsize)+iirht, hash);
}
}
if (nrhtval != -1) nrht[ic] = nrhtval;
}
if (nbotval == -1) {
// Taking care of boundary cells
// Force each boundary cell to point to itself on its boundary direction
if (jjcur < 1*IPOW2(levmx) -jminsize) nbotval = read_hash(jjcur*(imaxsize-iminsize)+iicur, hash);
// Boundary cells next to corner boundary need special checks
if (jjcur == 1*IPOW2(levmx)-jminsize && (iicur < 1*IPOW2(levmx)-iminsize || iicur >= imax*IPOW2(levmx)-iminsize ) ) nbotval = read_hash(jjcur*(imaxsize-iminsize)+iicur, hash);
// need to check for finer neighbor first
// Right and top neighbor don't change for finer, so drop through to same size
// Left and bottom need to be half of same size index for finer test
if (lev != levmx) {
int jjbotfiner = jjcur-(jjcur-jjbot)/2;
if (nbotval == -1 && jjbotfiner >= 0) nbotval = read_hash(jjbotfiner*(imaxsize-iminsize)+iicur, hash);
}
// same size neighbor
if (nbotval == -1 && jjbot >=0) nbotval = read_hash(jjbot*(imaxsize-iminsize)+iicur, hash);
// Now we need to take care of special case where bottom and left boundary need adjustment since
// expected cell doesn't exist on these boundaries if it is finer than current cell
if (iicur < 1*IPOW2(levmx) && lev != levmx) {
if (nbotval == -1) {
int iirhtfiner = (iicur+iirht)/2;
int jjbotfiner = jjcur-(jjcur-jjbot)/2;
if (jjbotfiner >= 0 && iirhtfiner < imaxsize-iminsize) nbotval = read_hash(jjbotfiner*(imaxsize-iminsize)+iirhtfiner, hash);
}
}
// coarser neighbor
if (lev != 0){
if (nbotval == -1) {
int jjbotcoarser = jjbot - (jjcur-jjbot);
int iibot = (ii/2)*2*levmult-iminsize;
if (jjbotcoarser >= 0 && iibot >= 0) nbotval = read_hash(jjbotcoarser*(imaxsize-iminsize)+iibot, hash);
}
}
if (nbotval != -1) nbot[ic] = nbotval;
}
if (ntopval == -1) {
// Taking care of boundary cells
// Force each boundary cell to point to itself on its boundary direction
if (jjcur > jmax*IPOW2(levmx)-1-jminsize) ntopval = read_hash(jjcur*(imaxsize-iminsize)+iicur, hash);
// Boundary cells next to corner boundary need special checks
if (jjtop == jmax*IPOW2(levmx)-jminsize && (iicur < 1*IPOW2(levmx)-iminsize || iicur >= imax*IPOW2(levmx)-iminsize ) ) ntopval = read_hash(jjcur*(imaxsize-iminsize)+iicur, hash);
// same size neighbor
if (ntopval == -1 && jjtop < jmaxsize-jminsize) ntopval = read_hash(jjtop*(imaxsize-iminsize)+iicur, hash);
if (iicur < 1*IPOW2(levmx)) {
if (ntopval == -1) {
int iirhtfiner = (iicur+iirht)/2;
if (jjtop < jmaxsize-jminsize && iirhtfiner < imaxsize-iminsize) ntopval = read_hash(jjtop*(imaxsize-iminsize)+iirhtfiner, hash);
}
}
// coarser neighbor
if (lev != 0){
if (ntopval == -1) {
int iitop = (ii/2)*2*levmult-iminsize;
if (jjtop < jmaxsize-jminsize && iitop < imaxsize-iminsize) ntopval = read_hash(jjtop*(imaxsize-iminsize)+iitop, hash);
}
}
if (ntopval != -1) ntop[ic] = ntopval;
}
//fprintf(fp,"%d: neighbors[%d] = %d %d %d %d\n",mype,ic,nlft[ic],nrht[ic],nbot[ic],ntop[ic]);
}
if (TIMING_LEVEL >= 2) {
#ifdef _OPENMP
#pragma omp master
#endif
cpu_timers[MESH_TIMER_FILL_NEIGH_GHOST] += cpu_timer_stop(tstart_lev2);
cpu_timer_start(&tstart_lev2);
}
if (DEBUG) {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
fprintf(fp,"After setting neighbors through ghost cells\n");
print_local();
#ifdef _OPENMP
} // end master region
#endif
}
/*
// Set neighbors to global cell numbers from hash
for (uint ic=0; ic<ncells; ic++){
ii = i[ic];
jj = j[ic];
lev = level[ic];
levmult = IPOW2(levmx-lev);
//fprintf(fp,"%d:Neighbors input for ic %d ii %d jj %d levmult %d lev %d\n",mype,ic, ii, jj, levmult,lev);
//fprintf(fp,"%d:Neighbors befor ic %d nlft %d nrht %d nbot %d ntop %d\n",mype,ic,nlft[ic],nrht[ic],nbot[ic],ntop[ic]);
if (nlft[ic] == -1) nlft[ic] = hash[( jj *levmult )-jminsize][(max( ii *levmult-1, 0 ))-iminsize];
if (celltype[ic] == BOTTOM_BOUNDARY && nlft[ic] == -1){
if (nlft[ic] == -1) nlft[ic] = hash[(jj+1)*levmult-jminsize][(min( (ii+1)*levmult, imaxcalc-1))-iminsize];
}
if (nrht[ic] == -1) nrht[ic] = hash[( jj *levmult )-jminsize][(min( (ii+1)*levmult, imaxcalc-1))-iminsize];
if (celltype[ic] == BOTTOM_BOUNDARY && nrht[ic] == -1){
if (nrht[ic] == -1) nrht[ic] = hash[(jj+1)*levmult-jminsize][(min( (ii+1)*levmult, imaxcalc-1))-iminsize];
//if (ic == 3 && mype == 0) printf("DEBUG line %d -- ic %d celltype %d nrht %d\n",__line__,ic,celltype[ic],nrht[ic]);
//printf("DEBUG line %d -- ic %d celltype %d nrht %d jj %d ii %d\n",__line__,ic,celltype[ic],nrht[ic],(jj+1)*levmult-jminsize,(min( (ii+1)*levmult, imaxcalc-1))-iminsize);
}
if (nbot[ic] == -1) nbot[ic] = hash[(max( jj *levmult-1, 0) )-jminsize][( ii *levmult )-iminsize];
if (celltype[ic] == LEFT_BOUNDARY && nbot[ic] == -1){
if (nbot[ic] == -1) nbot[ic] = hash[(max( jj *levmult-1, 0) )-jminsize][( ii *levmult+1 )-iminsize];
}
if (ntop[ic] == -1) ntop[ic] = hash[(min( (jj+1)*levmult, jmaxcalc-1))-jminsize][( ii *levmult )-iminsize];
if (celltype[ic] == LEFT_BOUNDARY && ntop[ic] == -1){
if (ntop[ic] == -1) ntop[ic] = hash[(min( (jj+1)*levmult, jmaxcalc-1))-jminsize][( ii *levmult+1 )-iminsize];
}
//fprintf(fp,"%d:Neighbors after ic %d nlft %d nrht %d nbot %d ntop %d\n",mype,ic,nlft[ic],nrht[ic],nbot[ic],ntop[ic]);
}
*/
if (TIMING_LEVEL >= 2) {
#ifdef _OPENMP
#pragma omp master
#endif
cpu_timers[MESH_TIMER_SET_CORNER_NEIGH] += cpu_timer_stop(tstart_lev2);
cpu_timer_start(&tstart_lev2);
}
if (DEBUG) {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
fprintf(fp,"After setting corner neighbors\n");
print_local();
#ifdef _OPENMP
} // end master region
#endif
}
// Adjusting neighbors to local indices
#ifdef _OPENMP
#pragma omp for
#endif
for (uint ic=0; ic<ncells_ghost; ic++){
//fprintf(fp,"%d: ic %d nlft %d noffset %d ncells %ld\n",mype,ic,nlft[ic],noffset,ncells);
if (nlft[ic] <= -(int)ncells && nlft[ic] > -(int)ncells_ghost){
nlft[ic] = abs(nlft[ic]);
} else if (nlft[ic] >= noffset && nlft[ic] < (int)(noffset+ncells)) {
nlft[ic] -= noffset;
}
if (nrht[ic] <= -(int)ncells && nrht[ic] > -(int)ncells_ghost){
nrht[ic] = abs(nrht[ic]);
} else if (nrht[ic] >= noffset && nrht[ic] < (int)(noffset+ncells)) {
nrht[ic] -= noffset;
}
if (nbot[ic] <= -(int)ncells && nbot[ic] > -(int)ncells_ghost){
nbot[ic] = abs(nbot[ic]);
} else if (nbot[ic] >= noffset && nbot[ic] < (int)(noffset+ncells)) {
nbot[ic] -= noffset;
}
if (ntop[ic] <= -(int)ncells && ntop[ic] > -(int)ncells_ghost){
ntop[ic] = abs(ntop[ic]);
} else if (ntop[ic] >= noffset && ntop[ic] < (int)(noffset+ncells)) {
ntop[ic] -= noffset;
}
}
if (DEBUG) {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
fprintf(fp,"After adjusting neighbors to local indices\n");
print_local();
#ifdef _OPENMP
} // end master region
#endif
}
if (TIMING_LEVEL >= 2) {
#ifdef _OPENMP
#pragma omp master
#endif
cpu_timers[MESH_TIMER_NEIGH_ADJUST] += cpu_timer_stop(tstart_lev2);
cpu_timer_start(&tstart_lev2);
}
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
offtile_ratio_local = (offtile_ratio_local*(double)offtile_local_count) + ((double)nghost / (double)ncells);
offtile_local_count++;
offtile_ratio_local /= offtile_local_count;
//printf("%d ncells size is %ld ncells_ghost size is %ld nghost %d\n",mype,ncells,ncells_ghost,nghost);
//fprintf(fp,"%d ncells_ghost size is %ld nghost %d\n",mype,ncells_ghost,nghost);
if (cell_handle) L7_Free(&cell_handle);
cell_handle=0;
if (DEBUG) {
fprintf(fp,"%d: SETUP ncells %ld noffset %d nghost %d\n",mype,ncells,noffset,nghost);
for (int ig = 0; ig<nghost; ig++){
fprintf(fp,"%d: indices needed ic %d index %d\n",mype,ig,indices_needed[ig]);
}
}
L7_Setup(0, noffset, ncells, &indices_needed[0], nghost, &cell_handle);
if (TIMING_LEVEL >= 2) cpu_timers[MESH_TIMER_SETUP_COMM] += cpu_timer_stop(tstart_lev2);
#ifdef _OPENMP
} // end master region
#endif
if (DEBUG) {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
print_local();
int jmaxglobal = (jmax+1)*IPOW2(levmx);
int imaxglobal = (imax+1)*IPOW2(levmx);
fprintf(fp,"\n HASH numbering\n");
for (int jj = jmaxglobal-1; jj>=0; jj--){
fprintf(fp,"%2d: %4d:",mype,jj);
if (jj >= jminsize && jj < jmaxsize) {
for (int ii = 0; ii<imaxglobal; ii++){
if (ii >= iminsize && ii < imaxsize) {
fprintf(fp,"%5d",read_hash((jj-jminsize)*(imaxsize-iminsize)+(ii-iminsize), hash) );
} else {
fprintf(fp," ");
}
}
}
fprintf(fp,"\n");
}
fprintf(fp,"%2d: ",mype);
for (int ii = 0; ii<imaxglobal; ii++){
fprintf(fp,"%4d:",ii);
}
fprintf(fp,"\n");
fprintf(fp,"\n nlft numbering\n");
for (int jj = jmaxglobal-1; jj>=0; jj--){
fprintf(fp,"%2d: %4d:",mype,jj);
if (jj >= jminsize && jj < jmaxsize) {
for (int ii = 0; ii<imaxglobal; ii++){
if (ii >= iminsize && ii < imaxsize) {
int hashval = read_hash((jj-jminsize)*(imaxsize-iminsize)+(ii-iminsize), hash) -noffset;
if ( (hashval >= 0 && hashval < (int)ncells) ) {
fprintf(fp,"%5d",nlft[hashval]);
} else {
fprintf(fp," ");
}
}
}
}
fprintf(fp,"\n");
}
fprintf(fp,"%2d: ",mype);
for (int ii = 0; ii<imaxglobal; ii++){
fprintf(fp,"%4d:",ii);
}
fprintf(fp,"\n");
fprintf(fp,"\n nrht numbering\n");
for (int jj = jmaxglobal-1; jj>=0; jj--){
fprintf(fp,"%2d: %4d:",mype,jj);
if (jj >= jminsize && jj < jmaxsize) {
for (int ii = 0; ii<imaxglobal; ii++){
if ( ii >= iminsize && ii < imaxsize ) {
int hashval = read_hash((jj-jminsize)*(imaxsize-iminsize)+(ii-iminsize), hash) -noffset;
if ( hashval >= 0 && hashval < (int)ncells ) {
fprintf(fp,"%5d",nrht[hashval]);
} else {
fprintf(fp," ");
}
}
}
}
fprintf(fp,"\n");
}
fprintf(fp,"%2d: ",mype);
for (int ii = 0; ii<imaxglobal; ii++){
fprintf(fp,"%4d:",ii);
}
fprintf(fp,"\n");
fprintf(fp,"\n nbot numbering\n");
for (int jj = jmaxglobal-1; jj>=0; jj--){
fprintf(fp,"%2d: %4d:",mype,jj);
if (jj >= jminsize && jj < jmaxsize) {
for (int ii = 0; ii<imaxglobal; ii++){
if ( ii >= iminsize && ii < imaxsize ) {
int hashval = read_hash((jj-jminsize)*(imaxsize-iminsize)+(ii-iminsize), hash) -noffset;
if ( hashval >= 0 && hashval < (int)ncells ) {
fprintf(fp,"%5d",nbot[hashval]);
} else {
fprintf(fp," ");
}
}
}
}
fprintf(fp,"\n");
}
fprintf(fp,"%2d: ",mype);
for (int ii = 0; ii<imaxglobal; ii++){
fprintf(fp,"%4d:",ii);
}
fprintf(fp,"\n");
fprintf(fp,"\n ntop numbering\n");
for (int jj = jmaxglobal-1; jj>=0; jj--){
fprintf(fp,"%2d: %4d:",mype,jj);
if (jj >= jminsize && jj < jmaxsize) {
for (int ii = 0; ii<imaxglobal; ii++){
if ( ii >= iminsize && ii < imaxsize ) {
int hashval = read_hash((jj-jminsize)*(imaxsize-iminsize)+(ii-iminsize), hash) -noffset;
if ( hashval >= 0 && hashval < (int)ncells ) {
fprintf(fp,"%5d",ntop[hashval]);
} else {
fprintf(fp," ");
}
}
}
}
fprintf(fp,"\n");
}
fprintf(fp,"%2d: ",mype);
for (int ii = 0; ii<imaxglobal; ii++){
fprintf(fp,"%4d:",ii);
}
fprintf(fp,"\n");
#ifdef _OPENMP
} // end master region
#endif
} // end DEBUG
if (DEBUG) {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
print_local();
for (uint ic=0; ic<ncells; ic++){
fprintf(fp,"%d: before update ic %d i %d j %d lev %d nlft %d nrht %d nbot %d ntop %d\n",
mype,ic,i[ic],j[ic],level[ic],nlft[ic],nrht[ic],nbot[ic],ntop[ic]);
}
int ig=0;
for (uint ic=ncells; ic<ncells_ghost; ic++, ig++){
fprintf(fp,"%d: after update ic %d off %d i %d j %d lev %d nlft %d nrht %d nbot %d ntop %d\n",
mype,ic,indices_needed[ig],i[ic],j[ic],level[ic],nlft[ic],nrht[ic],nbot[ic],ntop[ic]);
}
#ifdef _OPENMP
} // end master region
#endif
} // end DEBUG
} // if numpe > 1
#endif
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
write_hash_collision_report();
read_hash_collision_report();
compact_hash_delete(hash);
#ifdef BOUNDS_CHECK
{
for (uint ic=0; ic<ncells; ic++){
int nl = nlft[ic];
if (nl<0 || nl>= (int)ncells_ghost) printf("%d: Warning at line %d cell %d nlft %d\n",mype,__LINE__,ic,nl);
if (level[nl] > level[ic]){
int ntl = ntop[nl];
if (ntl<0 || ntl>= (int)ncells_ghost) printf("%d: Warning at line %d cell %d global %d nlft %d ntop of nlft %d\n",mype,__LINE__,ic,ic+noffset,nl,ntl);
}
int nr = nrht[ic];
if (nr<0 || nr>= (int)ncells_ghost) printf("%d: Warning at line %d cell %d nrht %d\n",mype,__LINE__,ic,nr);
if (level[nr] > level[ic]){
int ntr = ntop[nr];
if (ntr<0 || ntr>= (int)ncells_ghost) printf("%d: Warning at line %d cell %d ntop of nrht %d\n",mype,__LINE__,ic,ntr);
}
int nb = nbot[ic];
if (nb<0 || nb>= (int)ncells_ghost) printf("%d: Warning at line %d cell %d nbot %d\n",mype,__LINE__,ic,nb);
if (level[nb] > level[ic]){
int nrb = nrht[nb];
if (nrb<0 || nrb>= (int)ncells_ghost) printf("%d: Warning at line %d cell %d nrht of nbot %d\n",mype,__LINE__,ic,nrb);
}
int nt = ntop[ic];
if (nt<0 || nt>= (int)ncells_ghost) printf("%d: Warning at line %d cell %d ntop %d\n",mype,__LINE__,ic,nt);
if (level[nt] > level[ic]){
int nrt = nrht[nt];
if (nrt<0 || nrt>= (int)ncells_ghost) printf("%d: Warning at line %d cell %d nrht of ntop %d\n",mype,__LINE__,ic,nrt);
}
}
}
#endif
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
} else if (calc_neighbor_type == KDTREE) {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
struct timeval tstart_lev2;
if (TIMING_LEVEL >= 2) cpu_timer_start(&tstart_lev2);
TBounds box;
vector<int> index_list(IPOW2(levmx*levmx) );
int num;
ibase = 0;
calc_spatial_coordinates(ibase);
kdtree_setup();
if (TIMING_LEVEL >= 2) {
cpu_timers[MESH_TIMER_KDTREE_SETUP] += cpu_timer_stop(tstart_lev2);
cpu_timer_start(&tstart_lev2);
}
for (uint ic=0; ic<ncells; ic++) {
//left
nlft[ic] = ic;
box.min.x = x[ic]-0.25*dx[ic];
box.max.x = x[ic]-0.25*dx[ic];
box.min.y = y[ic]+0.25*dy[ic];
box.max.y = y[ic]+0.25*dy[ic];
KDTree_QueryBoxIntersect(&tree, &num, &(index_list[0]), &box);
if (num == 1) nlft[ic]=index_list[0];
//right
nrht[ic] = ic;
box.min.x = x[ic]+1.25*dx[ic];
box.max.x = x[ic]+1.25*dx[ic];
box.min.y = y[ic]+0.25*dy[ic];
box.max.y = y[ic]+0.25*dy[ic];
KDTree_QueryBoxIntersect(&tree, &num, &(index_list[0]), &box);
if (num == 1) nrht[ic]=index_list[0];
//bot
nbot[ic] = ic;
box.min.x = x[ic]+0.25*dx[ic];
box.max.x = x[ic]+0.25*dx[ic];
box.min.y = y[ic]-0.25*dy[ic];
box.max.y = y[ic]-0.25*dy[ic];
KDTree_QueryBoxIntersect(&tree, &num, &(index_list[0]), &box);
if (num == 1) nbot[ic]=index_list[0];
//top
ntop[ic] = ic;
box.min.x = x[ic]+0.25*dx[ic];
box.max.x = x[ic]+0.25*dx[ic];
box.min.y = y[ic]+1.25*dy[ic];
box.max.y = y[ic]+1.25*dy[ic];
KDTree_QueryBoxIntersect(&tree, &num, &(index_list[0]), &box);
if (num == 1) ntop[ic]=index_list[0];
} // End main loop over cells.
KDTree_Destroy(&tree);
if (TIMING_LEVEL >= 2) cpu_timers[MESH_TIMER_KDTREE_QUERY] += cpu_timer_stop(tstart_lev2);
#ifdef _OPENMP
}
#pragma omp barrier
#endif
} // calc_neighbor_type
}
#ifdef _OPENMP
#pragma omp master
#endif
cpu_timers[MESH_TIMER_CALC_NEIGHBORS] += cpu_timer_stop(tstart_cpu);
}