| /* |
| * Copyright (c) 2011-2012, Los Alamos National Security, LLC. |
| * All rights Reserved. |
| * |
| * Copyright 2011-2012. Los Alamos National Security, LLC. This software was produced |
| * under U.S. Government contract DE-AC52-06NA25396 for Los Alamos National |
| * Laboratory (LANL), which is operated by Los Alamos National Security, LLC |
| * for the U.S. Department of Energy. The U.S. Government has rights to use, |
| * reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR LOS |
| * ALAMOS NATIONAL SECURITY, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR |
| * ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is modified |
| * to produce derivative works, such modified software should be clearly marked, |
| * so as not to confuse it with the version available from LANL. |
| * |
| * Additionally, redistribution and use in source and binary forms, with or without |
| * modification, are permitted provided that the following conditions are met: |
| * * Redistributions of source code must retain the above copyright |
| * notice, this list of conditions and the following disclaimer. |
| * * Redistributions in binary form must reproduce the above copyright |
| * notice, this list of conditions and the following disclaimer in the |
| * documentation and/or other materials provided with the distribution. |
| * * Neither the name of the Los Alamos National Security, LLC, Los Alamos |
| * National Laboratory, LANL, the U.S. Government, nor the names of its |
| * contributors may be used to endorse or promote products derived from |
| * this software without specific prior written permission. |
| * |
| * THIS SOFTWARE IS PROVIDED BY THE LOS ALAMOS NATIONAL SECURITY, LLC AND |
| * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT |
| * NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
| * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL LOS ALAMOS NATIONAL |
| * SECURITY, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, |
| * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
| * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; |
| * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, |
| * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) |
| * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
| * POSSIBILITY OF SUCH DAMAGE. |
| * |
| * CLAMR -- LA-CC-11-094 |
| * This research code is being developed as part of the |
| * 2011 X Division Summer Workshop for the express purpose |
| * of a collaborative code for development of ideas in |
| * the implementation of AMR codes for Exascale platforms |
| * |
| * AMR implementation of the Wave code previously developed |
| * as a demonstration code for regular grids on Exascale platforms |
| * as part of the Supercomputing Challenge and Los Alamos |
| * National Laboratory |
| * |
| * Authors: Bob Robey XCP-2 brobey@lanl.gov |
| * Neal Davis davis68@lanl.gov, davis68@illinois.edu |
| * David Nicholaeff dnic@lanl.gov, mtrxknight@aol.com |
| * Dennis Trujillo dptrujillo@lanl.gov, dptru10@gmail.com |
| * |
| */ |
| #ifdef HAVE_MPI |
| #include "mpi.h" |
| #endif |
| |
| #include <sys/time.h> |
| #include <stdio.h> |
| #include <stdlib.h> |
| #include <math.h> |
| #include <list> |
| #include <algorithm> |
| #include "partition.h" |
| #include "KDTree.h" |
| #include "mesh.h" |
| #ifdef HAVE_MPI |
| #include "s7.h" |
| #endif |
| #include "zorder.h" |
| #include "timer.h" |
| #include "hsfc.h" |
| |
| #ifndef DEBUG |
| #define DEBUG 0 |
| #endif |
| |
| typedef unsigned int uint; |
| |
| int measure_type; |
| int meas_count = 0; |
| double meas_sum_average = 0.0; |
| |
| extern bool localStencil; |
| extern int initial_order; |
| extern int cycle_reorder; |
| |
| void Mesh::partition_measure(void) |
| { |
| if (measure_type != NO_PARTITION_MEASURE){ |
| |
| int ntX = TILE_SIZE; |
| static double offtile_ratio = 0.0; |
| |
| uint num_groups = (ncells + TILE_SIZE - 1)/TILE_SIZE; |
| |
| if (measure_type == WITH_DUPLICATES) { |
| int i = 0; |
| #ifdef _OPENMP |
| #pragma omp for reduction(+:offtile_ratio) |
| #endif |
| for (uint group_id=0; group_id < num_groups; group_id ++){ |
| |
| int start_idx = group_id * ntX; |
| int end_idx = (group_id + 1) * ntX; |
| |
| int offtile =0; |
| for (uint ic = 0; ic < TILE_SIZE; ic++, i++){ |
| |
| if (i >= ncells) continue; |
| //taken from wave_kern_calc.cl 'setup tile' kernel |
| if (nlft[i] < start_idx || nlft[i] >= end_idx) offtile++; |
| if (level[nlft[i]] > level[i] && |
| (ntop[nlft[i]] < start_idx || ntop[nlft[i]] >= end_idx) ) offtile++; |
| if (nrht[i] < start_idx || nrht[i] >= end_idx) offtile++; |
| if (level[nrht[i]] > level[i] && |
| (ntop[nrht[i]] < start_idx || ntop[nrht[i]] >= end_idx) ) offtile++; |
| if (nbot[i] < start_idx || nbot[i] >= end_idx) offtile++; |
| if (level[nbot[i]] > level[i] && |
| (nrht[nbot[i]] < start_idx || nrht[nbot[i]] >= end_idx) ) offtile++; |
| if (ntop[i] < start_idx || ntop[i] >= end_idx) offtile++; |
| if (level[ntop[i]] > level[i] && |
| (nrht[ntop[i]] < start_idx || nrht[ntop[i]] >= end_idx) ) offtile++; |
| } |
| offtile_ratio += (double)offtile/(double)(TILE_SIZE); |
| //printf("DEBUG Ratio of surface area to volume is equal to %d / %d ratio is %lf\n", offtile, TILE_SIZE, (double)offtile/(double)TILE_SIZE); |
| } |
| } else if (measure_type == WITHOUT_DUPLICATES) { |
| int i = 0; |
| #ifdef _OPENMP |
| #pragma omp for reduction(+:offtile_ratio) |
| #endif |
| for (uint group_id=0; group_id < num_groups; group_id ++){ |
| list<int> offtile_list; |
| |
| int start_idx = group_id * ntX; |
| int end_idx = (group_id + 1) * ntX; |
| |
| for (uint ic = 0; ic < TILE_SIZE; ic++, i++){ |
| |
| if (i >= ncells) continue; |
| |
| if (nlft[i] < start_idx || nlft[i] >= end_idx) offtile_list.push_back(nlft[i]); |
| if (level[nlft[i]] > level[i] && |
| (ntop[nlft[i]] < start_idx || ntop[nlft[i]] >= end_idx) ) offtile_list.push_back(ntop[nlft[i]]); |
| if (nrht[i] < start_idx || nrht[i] >= end_idx) offtile_list.push_back(nrht[i]); |
| if (level[nrht[i]] > level[i] && |
| (ntop[nrht[i]] < start_idx || ntop[nrht[i]] >= end_idx) ) offtile_list.push_back(ntop[nrht[i]]); |
| if (nbot[i] < start_idx || nbot[i] >= end_idx) offtile_list.push_back(nbot[i]); |
| if (level[nbot[i]] > level[i] && |
| (nrht[nbot[i]] < start_idx || nrht[nbot[i]] >= end_idx) ) offtile_list.push_back(nrht[nbot[i]]); |
| if (ntop[i] < start_idx || ntop[i] >= end_idx) offtile_list.push_back(ntop[i]); |
| if (level[ntop[i]] > level[i] && |
| (nrht[ntop[i]] < start_idx || nrht[ntop[i]] >= end_idx) ) offtile_list.push_back(nrht[ntop[i]]); |
| } |
| offtile_list.sort(); |
| offtile_list.unique(); |
| |
| offtile_ratio += (double)offtile_list.size()/(double)(TILE_SIZE); |
| //printf("DEBUG Ratio of surface area to volume is equal to %d / %d ratio is %lf\n", offtile, TILE_SIZE, (double)offtile/(double)TILE_SIZE); |
| } |
| } else if (measure_type == CVALUE) { |
| int i = 0; |
| #ifdef _OPENMP |
| #pragma omp for reduction(+:offtile_ratio) |
| #endif |
| for (uint group_id=0; group_id < num_groups; group_id ++){ |
| list<int> offtile_list; |
| |
| int start_idx = group_id * ntX; |
| int end_idx = (group_id + 1) * ntX; |
| |
| for (uint ic = 0; ic < TILE_SIZE; ic++, i++){ |
| |
| if (i >= ncells) continue; |
| |
| if (nlft[i] < start_idx || nlft[i] >= end_idx) offtile_list.push_back(nlft[i]); |
| if (level[nlft[i]] > level[i] && |
| (ntop[nlft[i]] < start_idx || ntop[nlft[i]] >= end_idx) ) offtile_list.push_back(ntop[nlft[i]]); |
| if (nrht[i] < start_idx || nrht[i] >= end_idx) offtile_list.push_back(nrht[i]); |
| if (level[nrht[i]] > level[i] && |
| (ntop[nrht[i]] < start_idx || ntop[nrht[i]] >= end_idx) ) offtile_list.push_back(ntop[nrht[i]]); |
| if (nbot[i] < start_idx || nbot[i] >= end_idx) offtile_list.push_back(nbot[i]); |
| if (level[nbot[i]] > level[i] && |
| (nrht[nbot[i]] < start_idx || nrht[nbot[i]] >= end_idx) ) offtile_list.push_back(nrht[nbot[i]]); |
| if (ntop[i] < start_idx || ntop[i] >= end_idx) offtile_list.push_back(ntop[i]); |
| if (level[ntop[i]] > level[i] && |
| (nrht[ntop[i]] < start_idx || nrht[ntop[i]] >= end_idx) ) offtile_list.push_back(nrht[ntop[i]]); |
| } |
| offtile_list.sort(); |
| offtile_list.unique(); |
| |
| offtile_ratio += (double)offtile_list.size()/(4*sqrt((double)(TILE_SIZE))); |
| //printf("DEBUG Ratio of surface area to volume is equal to %d / %d ratio is %lf\n", offtile, TILE_SIZE, (double)offtile/(double)TILE_SIZE); |
| } |
| } else if (measure_type == CSTARVALUE) { |
| int i = 0; |
| #ifdef _OPENMP |
| #pragma omp for reduction(+:offtile_ratio) |
| #endif |
| for (uint group_id=0; group_id < num_groups; group_id ++){ |
| list<int> offtile_list; |
| list<int> offtile_cache_lines; // Assumes memory is aligned |
| int cache_line_size = 4; // Some could be 8, or more? |
| |
| int start_idx = group_id * ntX; |
| int end_idx = (group_id + 1) * ntX; |
| |
| for (uint ic = 0; ic < TILE_SIZE; ic++, i++){ |
| |
| if (i >= ncells) continue; |
| |
| if (nlft[i] < start_idx || nlft[i] >= end_idx) { |
| offtile_list.push_back(nlft[i]); |
| offtile_cache_lines.push_back(nlft[i]/cache_line_size); |
| } |
| |
| if (level[nlft[i]] > level[i] && (ntop[nlft[i]] < start_idx || ntop[nlft[i]] >= end_idx) ) { |
| offtile_list.push_back(ntop[nlft[i]]); |
| offtile_cache_lines.push_back(ntop[nlft[i]]/cache_line_size); |
| } |
| if (nrht[i] < start_idx || nrht[i] >= end_idx) { |
| offtile_list.push_back(nrht[i]); |
| offtile_cache_lines.push_back(nrht[i]/cache_line_size); |
| } |
| if (level[nrht[i]] > level[i] && (ntop[nrht[i]] < start_idx || ntop[nrht[i]] >= end_idx) ) { |
| offtile_list.push_back(ntop[nrht[i]]); |
| offtile_cache_lines.push_back(ntop[nrht[i]]/cache_line_size); |
| } |
| if (nbot[i] < start_idx || nbot[i] >= end_idx) { |
| offtile_list.push_back(nbot[i]); |
| offtile_cache_lines.push_back(nbot[i]/cache_line_size); |
| } |
| if (level[nbot[i]] > level[i] && (nrht[nbot[i]] < start_idx || nrht[nbot[i]] >= end_idx) ) { |
| offtile_list.push_back(nrht[nbot[i]]); |
| offtile_cache_lines.push_back(nrht[nbot[i]]/cache_line_size); |
| } |
| if (ntop[i] < start_idx || ntop[i] >= end_idx) { |
| offtile_list.push_back(ntop[i]); |
| offtile_cache_lines.push_back(ntop[i]/cache_line_size); |
| } |
| if (level[ntop[i]] > level[i] && (nrht[ntop[i]] < start_idx || nrht[ntop[i]] >= end_idx) ) { |
| offtile_list.push_back(nrht[ntop[i]]); |
| offtile_cache_lines.push_back(nrht[ntop[i]]/cache_line_size); |
| } |
| } |
| offtile_list.sort(); |
| offtile_list.unique(); |
| offtile_cache_lines.sort(); |
| offtile_cache_lines.unique(); |
| |
| double s_ngeom = (double)(offtile_list.size()); |
| double q_ngeom = (double)(offtile_cache_lines.size()); |
| double ngeom = (double)(TILE_SIZE); |
| double cover = (double)(cache_line_size); |
| // offtile_ratio += (s_ngeom * q_ngeom) / (4*sqrt(ngeom)*2*(1+(ngeom+cache_line_size-1)/cache_line_size)); |
| // offtile_ratio += (q_ngeom) / (2*sqrt(ngeom)+2*((sqrt(ngeom)+cover-1)/cover)); |
| // offtile_ratio += (q_ngeom) / ( (8*sqrt(ngeom)+cover-1)/cover ); |
| ngeom = sqrt(ngeom); |
| offtile_ratio += (s_ngeom*q_ngeom*cover) / ( 4 * ngeom * (8*ngeom+cover-1) ); |
| |
| //printf("DEBUG Ratio of surface area to volume is equal to %d / %d ratio is %lf\n", offtile, TILE_SIZE, (double)offtile/(double)TILE_SIZE); |
| } |
| } |
| |
| // printf("DEBUG Ratio of surface area to volume is equal to %d / %d \n", offtile, ontile); |
| |
| #ifdef _OPENMP |
| #pragma omp master |
| { |
| #endif |
| meas_count ++; |
| meas_sum_average += offtile_ratio/(double)num_groups; |
| // printf("DEBUG %d icount %d sum_average %lf\n",__LINE__,icount, sum_average); |
| #ifdef _OPENMP |
| } |
| #endif |
| } // if PARTITION TYPE |
| } |
| |
| void Mesh::print_partition_measure() |
| { |
| if (meas_count != 0) { |
| if (measure_type == NO_PARTITION_MEASURE) { |
| if (mype == 0) printf("No Partition Measure\n"); |
| } else if (measure_type == WITH_DUPLICATES) { |
| parallel_output("Average surface area to volume ratio ", meas_sum_average/(double)meas_count, 0, "with duplicates"); |
| } else if (measure_type == WITHOUT_DUPLICATES) { |
| parallel_output("Average surface area to volume ratio ", meas_sum_average/(double)meas_count, 0, "without duplicates"); |
| } else if (measure_type == CVALUE) { |
| parallel_output("Partition Quality Avg C value ", meas_sum_average/(double)meas_count, 0, ""); |
| } else if (measure_type == CSTARVALUE){ |
| parallel_output("Partition Quality Avg C* value ", meas_sum_average/(double)meas_count, 0, ""); |
| } |
| } |
| |
| if (numpe > 1){ |
| parallel_output("The MPI surface area to volume ratio ", offtile_ratio_local, 0, "without duplicates"); |
| } |
| } |
| |
| void Mesh::print_partition_type() |
| { |
| if (mype == 0) { |
| if (initial_order == ORIGINAL_ORDER) { |
| printf("Initial order is naive."); |
| } else if (initial_order == HILBERT_SORT) { |
| printf("Initial order is Hilbert sort."); |
| } else if (initial_order == HILBERT_PARTITION) { |
| printf("Initial order is Hilbert partitionr."); |
| } else if (initial_order == ZORDER) { |
| printf("Initial order is Z order."); |
| } |
| |
| if (cycle_reorder == ORIGINAL_ORDER) { |
| printf(" No cycle reorder."); |
| } else if (cycle_reorder == HILBERT_SORT) { |
| printf(" Cycle reorder is Hilbert sort."); |
| } else if (cycle_reorder == HILBERT_PARTITION) { |
| printf(" Cycle reorder is Hilbert partition."); |
| } else if (cycle_reorder == ZORDER) { |
| printf(" Cycle reorder is Z order."); |
| } |
| |
| if (localStencil) { |
| printf(" Local Stencil is on.\n"); |
| } else { |
| printf("\n"); |
| } |
| } |
| |
| } |
| void Mesh::partition_cells( |
| int numpe, // |
| vector<int> &z_order, // Resulting index ordering. |
| enum partition_method method) // Assigned partitioning method. |
| { |
| int *info; // |
| double iscale, // |
| jscale; // |
| int imax, // Maximum x-index. |
| jmax; // Maximum y-index. |
| vector<int> z_index; // Ordered curve from hsfc. |
| vector<int> i_scaled; // x-indices normalized to a scale of [0, 1] for hsfc. |
| vector<int> j_scaled; // y-indices normalized to a scale of [0, 1] for hsfc. |
| vector<double> iunit; // |
| vector<double> junit; // |
| |
| struct timeval tstart_cpu; |
| cpu_timer_start(&tstart_cpu); |
| |
| // Initialize ordered curve index. |
| z_index.resize(ncells, 0); |
| //z_order.resize(ncells, 0); |
| |
| if (parallel) { |
| #ifdef HAVE_MPI |
| nsizes.resize(numpe); |
| ndispl.resize(numpe); |
| MPI_Allgather(&ncells, 1, MPI_INT, &nsizes[0], 1, MPI_INT, MPI_COMM_WORLD); |
| ndispl[0]=0; |
| for (int ip=1; ip<numpe; ip++){ |
| ndispl[ip] = ndispl[ip-1] + nsizes[ip-1]; |
| } |
| noffset=0; |
| for (int ip=0; ip<mype; ip++){ |
| noffset += nsizes[ip]; |
| } |
| #endif |
| } else { |
| // Adjust the number of required work items to the number of cells. |
| proc.resize(ncells); |
| // Decompose the domain equitably. |
| calc_distribution(numpe); |
| noffset = 0; |
| } |
| |
| |
| // Partition cells according to one of several possible orderings. |
| int have_spatial_variables=0; |
| switch (method) |
| { case ORIGINAL_ORDER: |
| // Set z_order to the current cell order. |
| for (uint ic = 0; ic < ncells; ++ic) |
| { z_order[ic] = ic; } |
| |
| cpu_timers[MESH_TIMER_PARTITION] += cpu_timer_stop(tstart_cpu); |
| |
| return; |
| break; |
| |
| case HILBERT_SORT: |
| // Resort the curve by Hilbert order. |
| have_spatial_variables = 1; |
| if (x.size() < ncells) { |
| calc_spatial_coordinates(0); |
| have_spatial_variables = 0; |
| } |
| calc_centerminmax(); |
| iunit.resize(ncells); |
| junit.resize(ncells); |
| |
| // Get the range of values in the x- and y-directions and make the scale square. |
| iscale = 1.0 / (xcentermax - xcentermin); |
| jscale = 1.0 / (ycentermax - ycentermin); |
| |
| // Scale the indices to a normalized [0, 1] range for hsfc. |
| for (uint ic = 0; ic < ncells; ++ic){ |
| iunit[ic] = (x[ic] + 0.5 * dx[ic] - xcentermin) * iscale; |
| junit[ic] = (y[ic] + 0.5 * dy[ic] - ycentermin) * jscale; |
| } |
| |
| if (have_spatial_variables == 0){ |
| x.clear(); |
| dx.clear(); |
| y.clear(); |
| dy.clear(); |
| } |
| |
| if (parallel){ |
| #ifdef HAVE_MPI |
| info = (int *)malloc(sizeof(int) * 3 * ncells_global); |
| vector<double>iunit_global(ncells_global); |
| vector<double>junit_global(ncells_global); |
| vector<int>z_order_global(ncells_global); |
| |
| MPI_Allgatherv(&iunit[0], ncells, MPI_DOUBLE, &iunit_global[0], &nsizes[0], &ndispl[0], MPI_DOUBLE, MPI_COMM_WORLD); |
| MPI_Allgatherv(&junit[0], ncells, MPI_DOUBLE, &junit_global[0], &nsizes[0], &ndispl[0], MPI_DOUBLE, MPI_COMM_WORLD); |
| // Sort the mesh into an ordered space-filling curve from hsfc. |
| hsfc2sort(ncells_global, &iunit_global[0], &junit_global[0], 0, info, 1); |
| |
| // Copy the cell order information from info into z_order. |
| for (uint ic = 0; ic < ncells_global; ++ic) |
| { z_order_global[ic] = info[ic]; } |
| free(info); |
| |
| // Order the mesh according to the calculated order (note that z_order is for both curves). |
| vector<int> int_global(ncells_global); |
| vector<int> int_global_new(ncells_global); |
| |
| // gather, reorder and scatter i |
| MPI_Allgatherv(&i[0], ncells, MPI_INT, &int_global[0], &nsizes[0], &ndispl[0], MPI_INT, MPI_COMM_WORLD); |
| for (int ic = 0; ic<(int)ncells_global; ic++){ |
| int_global_new[ic] = int_global[z_order_global[ic]]; |
| } |
| MPI_Scatterv(&int_global_new[0], &nsizes[0], &ndispl[0], MPI_INT, &i[0], ncells, MPI_INT, 0, MPI_COMM_WORLD); |
| |
| // gather, reorder and scatter j |
| MPI_Allgatherv(&j[0], ncells, MPI_INT, &int_global[0], &nsizes[0], &ndispl[0], MPI_INT, MPI_COMM_WORLD); |
| for (int ic = 0; ic<(int)ncells_global; ic++){ |
| int_global_new[ic] = int_global[z_order_global[ic]]; |
| } |
| MPI_Scatterv(&int_global_new[0], &nsizes[0], &ndispl[0], MPI_INT, &j[0], ncells, MPI_INT, 0, MPI_COMM_WORLD); |
| |
| // gather, reorder and scatter level |
| MPI_Allgatherv(&level[0], ncells, MPI_INT, &int_global[0], &nsizes[0], &ndispl[0], MPI_INT, MPI_COMM_WORLD); |
| for (int ic = 0; ic<(int)ncells_global; ic++){ |
| int_global_new[ic] = int_global[z_order_global[ic]]; |
| } |
| MPI_Scatterv(&int_global_new[0], &nsizes[0], &ndispl[0], MPI_INT, &level[0], ncells, MPI_INT, 0, MPI_COMM_WORLD); |
| |
| // It is faster just to recalculate these variables instead of communicating them |
| if (mesh_memory.get_memory_size(celltype) >= ncells) { |
| calc_celltype(mesh_memory.get_memory_size(celltype)); |
| } |
| |
| if (have_spatial_variables) { |
| calc_spatial_coordinates(0); |
| } |
| |
| if (mesh_memory.get_memory_size(nlft) >= ncells) { |
| vector<int> inv_z_order(ncells_global); |
| for (int ic = 0; ic<(int)ncells_global; ic++){ |
| inv_z_order[z_order_global[ic]] = ic; |
| } |
| |
| MPI_Allgatherv(&nlft[0], ncells, MPI_INT, &int_global[0], &nsizes[0], &ndispl[0], MPI_INT, MPI_COMM_WORLD); |
| for (int ic = 0; ic<(int)ncells_global; ic++){ |
| int_global_new[ic] = inv_z_order[int_global[z_order_global[ic]]]; |
| } |
| MPI_Scatterv(&int_global_new[0], &nsizes[0], &ndispl[0], MPI_INT, &nlft[0], ncells, MPI_INT, 0, MPI_COMM_WORLD); |
| |
| MPI_Allgatherv(&nrht[0], ncells, MPI_INT, &int_global[0], &nsizes[0], &ndispl[0], MPI_INT, MPI_COMM_WORLD); |
| for (int ic = 0; ic<(int)ncells_global; ic++){ |
| int_global_new[ic] = inv_z_order[int_global[z_order_global[ic]]]; |
| } |
| MPI_Scatterv(&int_global_new[0], &nsizes[0], &ndispl[0], MPI_INT, &nrht[0], ncells, MPI_INT, 0, MPI_COMM_WORLD); |
| |
| MPI_Allgatherv(&nbot[0], ncells, MPI_INT, &int_global[0], &nsizes[0], &ndispl[0], MPI_INT, MPI_COMM_WORLD); |
| for (int ic = 0; ic<(int)ncells_global; ic++){ |
| int_global_new[ic] = inv_z_order[int_global[z_order_global[ic]]]; |
| } |
| MPI_Scatterv(&int_global_new[0], &nsizes[0], &ndispl[0], MPI_INT, &nbot[0], ncells, MPI_INT, 0, MPI_COMM_WORLD); |
| |
| MPI_Allgatherv(&ntop[0], ncells, MPI_INT, &int_global[0], &nsizes[0], &ndispl[0], MPI_INT, MPI_COMM_WORLD); |
| for (int ic = 0; ic<(int)ncells_global; ic++){ |
| int_global_new[ic] = inv_z_order[int_global[z_order_global[ic]]]; |
| } |
| MPI_Scatterv(&int_global_new[0], &nsizes[0], &ndispl[0], MPI_INT, &ntop[0], ncells, MPI_INT, 0, MPI_COMM_WORLD); |
| } |
| |
| MPI_Scatterv(&z_order_global[0], &nsizes[0], &ndispl[0], MPI_INT, &z_order[0], ncells, MPI_INT, 0, MPI_COMM_WORLD); |
| #endif |
| } else { |
| info = (int *)malloc(sizeof(int) * 3 * ncells); |
| |
| // Sort the mesh into an ordered space-filling curve from hsfc. |
| hsfc2sort(ncells, &iunit[0], &junit[0], 0, info, 1); |
| |
| // Copy the cell order information from info into z_order. |
| for (uint ic = 0; ic < ncells; ++ic) |
| { z_order[ic] = info[ic]; } |
| free(info); |
| |
| // Order the mesh according to the calculated order (note that z_order is for both curves). |
| vector<int> int_local(ncells); |
| |
| mesh_memory.set_memory_attribute(nlft, 0x100); |
| mesh_memory.set_memory_attribute(nrht, 0x100); |
| mesh_memory.set_memory_attribute(nbot, 0x100); |
| mesh_memory.set_memory_attribute(ntop, 0x100); |
| |
| mesh_memory.memory_reorder_all(&z_order[0]); |
| memory_reset_ptrs(); |
| |
| if (x.size() >= ncells) { |
| vector<spatial_t> real_local(ncells); |
| |
| for (int ic = 0; ic<(int)ncells; ic++){ |
| real_local[ic] = x[ic]; |
| } |
| for (int ic = 0; ic<(int)ncells; ic++){ |
| x[ic] = real_local[z_order[ic]]; |
| } |
| |
| for (int ic = 0; ic<(int)ncells; ic++){ |
| real_local[ic] = dx[ic]; |
| } |
| for (int ic = 0; ic<(int)ncells; ic++){ |
| dx[ic] = real_local[z_order[ic]]; |
| } |
| |
| for (int ic = 0; ic<(int)ncells; ic++){ |
| real_local[ic] = y[ic]; |
| } |
| for (int ic = 0; ic<(int)ncells; ic++){ |
| y[ic] = real_local[z_order[ic]]; |
| } |
| |
| for (int ic = 0; ic<(int)ncells; ic++){ |
| real_local[ic] = dy[ic]; |
| } |
| for (int ic = 0; ic<(int)ncells; ic++){ |
| dy[ic] = real_local[z_order[ic]]; |
| } |
| } |
| |
| } |
| |
| break; |
| |
| case ZORDER: |
| // Resort the curve by z-order. |
| if (parallel) { |
| #ifdef HAVE_MPI |
| vector<int>i_global(ncells_global); |
| vector<int>j_global(ncells_global); |
| vector<int>level_global(ncells_global); |
| vector<int>z_index_global(ncells_global); |
| vector<int>z_order_global(ncells_global); |
| MPI_Allgatherv(&i[0], ncells, MPI_REAL, &i_global[0], &nsizes[0], &ndispl[0], MPI_REAL, MPI_COMM_WORLD); |
| MPI_Allgatherv(&j[0], ncells, MPI_REAL, &j_global[0], &nsizes[0], &ndispl[0], MPI_REAL, MPI_COMM_WORLD); |
| MPI_Allgatherv(&level[0], ncells, MPI_REAL, &level_global[0], &nsizes[0], &ndispl[0], MPI_REAL, MPI_COMM_WORLD); |
| |
| i_scaled.resize(ncells_global); |
| j_scaled.resize(ncells_global); |
| |
| // |
| imax = 0; |
| jmax = 0; |
| for (uint ic = 0; ic < ncells_global; ++ic) |
| { if (i_global[ic] > imax) imax = i_global[ic]; |
| if (j_global[ic] > jmax) jmax = j_global[ic]; } |
| |
| // |
| iscale = 16.0 / (double)imax; |
| jscale = 16.0 / (double)jmax; |
| |
| // |
| for (uint ic = 0; ic < ncells_global; ++ic) |
| { i_scaled[ic]=(int) ( (double)i_global[ic]*iscale); |
| j_scaled[ic]=(int) ( (double)j_global[ic]*jscale); } |
| |
| // |
| calc_zorder(ncells_global, &i_scaled[0], &j_scaled[0], &level_global[0], levmx, ibase, &z_index_global[0], &z_order_global[0]); |
| |
| // Order the mesh according to the calculated order (note that z_order is for both curves). |
| vector<int> int_global(ncells_global); |
| vector<int> int_global_new(ncells_global); |
| |
| // gather, reorder and scatter i |
| MPI_Allgatherv(&i[0], ncells, MPI_INT, &int_global[0], &nsizes[0], &ndispl[0], MPI_INT, MPI_COMM_WORLD); |
| for (int ic = 0; ic<(int)ncells_global; ic++){ |
| int_global_new[ic] = int_global[z_order_global[ic]]; |
| } |
| MPI_Scatterv(&int_global_new[0], &nsizes[0], &ndispl[0], MPI_INT, &i[0], ncells, MPI_INT, 0, MPI_COMM_WORLD); |
| |
| // gather, reorder and scatter j |
| MPI_Allgatherv(&j[0], ncells, MPI_INT, &int_global[0], &nsizes[0], &ndispl[0], MPI_INT, MPI_COMM_WORLD); |
| for (int ic = 0; ic<(int)ncells_global; ic++){ |
| int_global_new[ic] = int_global[z_order_global[ic]]; |
| } |
| MPI_Scatterv(&int_global_new[0], &nsizes[0], &ndispl[0], MPI_INT, &j[0], ncells, MPI_INT, 0, MPI_COMM_WORLD); |
| |
| // gather, reorder and scatter level |
| MPI_Allgatherv(&level[0], ncells, MPI_INT, &int_global[0], &nsizes[0], &ndispl[0], MPI_INT, MPI_COMM_WORLD); |
| for (int ic = 0; ic<(int)ncells_global; ic++){ |
| int_global_new[ic] = int_global[z_order_global[ic]]; |
| } |
| MPI_Scatterv(&int_global_new[0], &nsizes[0], &ndispl[0], MPI_INT, &level[0], ncells, MPI_INT, 0, MPI_COMM_WORLD); |
| |
| // It is faster just to recalculate these variables instead of communicating them |
| if (mesh_memory.get_memory_size(celltype) >= ncells) { |
| calc_celltype(mesh_memory.get_memory_size(celltype)); |
| } |
| |
| if (x.size() >= ncells) { |
| calc_spatial_coordinates(0); |
| } |
| |
| if (mesh_memory.get_memory_size(nlft) >= ncells) { |
| vector<int> inv_z_order(ncells_global); |
| for (int ic = 0; ic<(int)ncells_global; ic++){ |
| inv_z_order[z_order_global[ic]] = ic; |
| } |
| |
| MPI_Allgatherv(&nlft[0], ncells, MPI_INT, &int_global[0], &nsizes[0], &ndispl[0], MPI_INT, MPI_COMM_WORLD); |
| for (int ic = 0; ic<(int)ncells_global; ic++){ |
| int_global_new[ic] = inv_z_order[int_global[z_order_global[ic]]]; |
| } |
| MPI_Scatterv(&int_global_new[0], &nsizes[0], &ndispl[0], MPI_INT, &nlft[0], ncells, MPI_INT, 0, MPI_COMM_WORLD); |
| |
| MPI_Allgatherv(&nrht[0], ncells, MPI_INT, &int_global[0], &nsizes[0], &ndispl[0], MPI_INT, MPI_COMM_WORLD); |
| for (int ic = 0; ic<(int)ncells_global; ic++){ |
| int_global_new[ic] = inv_z_order[int_global[z_order_global[ic]]]; |
| } |
| MPI_Scatterv(&int_global_new[0], &nsizes[0], &ndispl[0], MPI_INT, &nrht[0], ncells, MPI_INT, 0, MPI_COMM_WORLD); |
| |
| MPI_Allgatherv(&nbot[0], ncells, MPI_INT, &int_global[0], &nsizes[0], &ndispl[0], MPI_INT, MPI_COMM_WORLD); |
| for (int ic = 0; ic<(int)ncells_global; ic++){ |
| int_global_new[ic] = inv_z_order[int_global[z_order_global[ic]]]; |
| } |
| MPI_Scatterv(&int_global_new[0], &nsizes[0], &ndispl[0], MPI_INT, &nbot[0], ncells, MPI_INT, 0, MPI_COMM_WORLD); |
| |
| MPI_Allgatherv(&ntop[0], ncells, MPI_INT, &int_global[0], &nsizes[0], &ndispl[0], MPI_INT, MPI_COMM_WORLD); |
| for (int ic = 0; ic<(int)ncells_global; ic++){ |
| int_global_new[ic] = inv_z_order[int_global[z_order_global[ic]]]; |
| } |
| MPI_Scatterv(&int_global_new[0], &nsizes[0], &ndispl[0], MPI_INT, &ntop[0], ncells, MPI_INT, 0, MPI_COMM_WORLD); |
| } |
| MPI_Scatterv(&z_order_global[0], &nsizes[0], &ndispl[0], MPI_REAL, &z_order[0], ncells, MPI_REAL, 0, MPI_COMM_WORLD); |
| #endif |
| } else { |
| i_scaled.resize(ncells); |
| j_scaled.resize(ncells); |
| |
| // |
| imax = 0; |
| jmax = 0; |
| for (uint ic = 0; ic < ncells; ++ic) |
| { if (i[ic] > imax) imax = i[ic]; |
| if (j[ic] > jmax) jmax = j[ic]; } |
| |
| // |
| iscale = 16.0 / (double)imax; |
| jscale = 16.0 / (double)jmax; |
| |
| // |
| for (uint ic = 0; ic < ncells; ++ic) |
| { i_scaled[ic]=(int) ( (double)i[ic]*iscale); |
| j_scaled[ic]=(int) ( (double)j[ic]*jscale); } |
| |
| // |
| calc_zorder(ncells, &i_scaled[0], &j_scaled[0], &level[0], levmx, ibase, &z_index[0], &z_order[0]); |
| |
| // Order the mesh according to the calculated order (note that z_order is for both curves). |
| vector<int> int_local(ncells); |
| |
| mesh_memory.set_memory_attribute(nlft, 0x100); |
| mesh_memory.set_memory_attribute(nrht, 0x100); |
| mesh_memory.set_memory_attribute(nbot, 0x100); |
| mesh_memory.set_memory_attribute(ntop, 0x100); |
| |
| mesh_memory.memory_reorder_all(&z_order[0]); |
| memory_reset_ptrs(); |
| |
| |
| if (x.size() >= ncells) { |
| vector<spatial_t> real_local(ncells); |
| |
| for (int ic = 0; ic<(int)ncells; ic++){ |
| real_local[ic] = x[ic]; |
| } |
| for (int ic = 0; ic<(int)ncells; ic++){ |
| x[ic] = real_local[z_order[ic]]; |
| } |
| |
| for (int ic = 0; ic<(int)ncells; ic++){ |
| real_local[ic] = dx[ic]; |
| } |
| for (int ic = 0; ic<(int)ncells; ic++){ |
| dx[ic] = real_local[z_order[ic]]; |
| } |
| |
| for (int ic = 0; ic<(int)ncells; ic++){ |
| real_local[ic] = y[ic]; |
| } |
| for (int ic = 0; ic<(int)ncells; ic++){ |
| y[ic] = real_local[z_order[ic]]; |
| } |
| |
| for (int ic = 0; ic<(int)ncells; ic++){ |
| real_local[ic] = dy[ic]; |
| } |
| for (int ic = 0; ic<(int)ncells; ic++){ |
| dy[ic] = real_local[z_order[ic]]; |
| } |
| } |
| |
| } |
| |
| break; |
| |
| default: |
| // Note that HILBERT_PARTITION is not currently supported due to redundancy with HILBERT_SORT. |
| break; |
| } |
| |
| |
| // Output ordered mesh information. |
| if (DEBUG) |
| { printf("orig index i j lev nlft nrht nbot ntop xlow xhigh ylow yhigh z index z order\n"); |
| for (uint ic=0; ic<ncells; ic++){ |
| printf(" %6d %4d %4d %4d %4d %4d %4d %4d ", index[ic], j[ic], i[ic], level[ic], nlft[ic], nrht[ic], nbot[ic], ntop[ic]); |
| printf(" %8.2lf %8.2lf %8.2lf %8.2lf", x[ic], x[ic]+dx[ic], y[ic], y[ic]+dy[ic]); |
| printf(" %6d %5d\n", z_index[ic], z_order[ic]); } } |
| |
| cpu_timers[MESH_TIMER_PARTITION] += cpu_timer_stop(tstart_cpu); |
| } |
| |
| // The distribution needs to be modified in order to spread out extra cells equitably among the work items. |
| void Mesh::calc_distribution(int numpe) |
| { |
| uint lsize = 0; // |
| for (int ip = 0; ip < numpe; ++ip) { |
| lsize += proc.size()/numpe; |
| if (ip < (int)proc.size()%numpe) lsize++; |
| for (int ic = 0; ic < (int)lsize; ic++) { |
| proc[ic] = ip; |
| } |
| } |
| } |
| |
