diff --git a/Makefile b/Makefile index 5550625e26..ddb5ab10e3 100644 --- a/Makefile +++ b/Makefile @@ -46,6 +46,7 @@ COMPFLAGS += -D OMPI_SKIP_MPICXX COMPFLAGS += -DPROFILE #Add -DNDEBUG to turn debugging off. If debugging is enabled performance will degrade significantly +#COMPFLAGS += -DDEBUG COMPFLAGS += -DNDEBUG # COMPFLAGS += -DIONOSPHERE_SORTED_SUMS # COMPFLAGS += -DDEBUG_SOLVERS diff --git a/grid.cpp b/grid.cpp index 14d7b349f0..04d69af269 100644 --- a/grid.cpp +++ b/grid.cpp @@ -115,24 +115,6 @@ void initializeGrids( MPI_Comm comm = MPI_COMM_WORLD; int neighborhood_size = VLASOV_STENCIL_WIDTH; - if (P::amrMaxSpatialRefLevel > 0) { - switch (VLASOV_STENCIL_WIDTH) { - case 1: - // Required cells will be included already - break; - case 2: - // looking from high to low refinement: stencil 2 will only give 1 cell, so need to add 1 - neighborhood_size = VLASOV_STENCIL_WIDTH+1; - break; - case 3: - // looking from high to low refinement: stencil 3 will only give 2 cells, so need to add 2 - // to reach surely into the third low-refinement neighbour - neighborhood_size = VLASOV_STENCIL_WIDTH+2; - break; - default: - std::cerr<<"Warning: unrecognized VLASOV_STENCIL_WIDTH in grid.cpp"< grid_length = {{P::xcells_ini, P::ycells_ini, P::zcells_ini}}; @@ -233,7 +215,7 @@ void initializeGrids( if (P::forceRefinement) { phiprof::Timer timer {"Restart refinement"}; - for (uint i = 0; i < P::amrMaxSpatialRefLevel; ++i) { + for (int i = 0; i < P::amrMaxSpatialRefLevel; ++i) { if (!adaptRefinement(mpiGrid, technicalGrid, sysBoundaries, project, i)) { cerr << "(MAIN) ERROR: Forcing refinement takes too much memory" << endl; exit(1); @@ -242,7 +224,7 @@ void initializeGrids( } } else if (P::refineOnRestart) { phiprof::Timer timer {"Restart refinement"}; - for (uint i = 0; i < P::amrMaxSpatialRefLevel; ++i) { + for (int i = 0; i < P::amrMaxSpatialRefLevel; ++i) { adaptRefinement(mpiGrid, technicalGrid, sysBoundaries, project); balanceLoad(mpiGrid, sysBoundaries); } @@ -1007,24 +989,6 @@ void initializeStencils(dccrg::Dccrg& mpi // In spatial AMR using DCCRG, the neighbors are considered relative to a given cell's size. // To get two coarse neighbors from a fine cell at interfaces, the stencil size needs to be increased by one. int addStencilDepth = 0; - if (P::amrMaxSpatialRefLevel > 0) { - switch (VLASOV_STENCIL_WIDTH) { - case 1: - // Required cells will be included already - break; - case 2: - // looking from high to low refinement: stencil 2 will only give 1 cell, so need to add 1 - addStencilDepth = 1; - break; - case 3: - // looking from high to low refinement: stencil 3 will only give 2 cells, so need to add 2 - // to reach surely into the third low-refinement neighbour - addStencilDepth = 2; - break; - default: - std::cerr<<"Warning: unrecognized VLASOV_STENCIL_WIDTH in grid.cpp"<& mpi } } } - /* Add extra face neighbors if required by AMR */ - for (int d = full_neighborhood_size+1; d <= full_neighborhood_size+addStencilDepth; d++) { - neighborhood.push_back({{ d, 0, 0}}); - neighborhood.push_back({{-d, 0, 0}}); - neighborhood.push_back({{0, d, 0}}); - neighborhood.push_back({{0,-d, 0}}); - neighborhood.push_back({{0, 0, d}}); - neighborhood.push_back({{0, 0,-d}}); - } /*all possible communication pairs*/ if( !mpiGrid.add_neighborhood(FULL_NEIGHBORHOOD_ID, neighborhood)){ std::cerr << "Failed to add neighborhood FULL_NEIGHBORHOOD_ID \n"; diff --git a/iowrite.cpp b/iowrite.cpp index 046c858ae9..1dcb267e2f 100644 --- a/iowrite.cpp +++ b/iowrite.cpp @@ -1104,7 +1104,7 @@ bool writeVelocitySpace(dccrg::Dccrg& mpi // Loop over AMR levels uint startindex=1; uint endindex=1; - for (uint AMR = 0; AMR <= P::amrMaxSpatialRefLevel; AMR++) { + for (int AMR = 0; AMR <= P::amrMaxSpatialRefLevel; AMR++) { int AMRm = 1u << AMR; uint cellsthislevel = (AMRm * P::xcells_ini) * (AMRm*P::ycells_ini) * (AMRm * P::zcells_ini); startindex = endindex; diff --git a/memoryallocation.h b/memoryallocation.h index 7f2974c11c..cb1544bdf4 100644 --- a/memoryallocation.h +++ b/memoryallocation.h @@ -24,6 +24,7 @@ #include #include +#include #include #include diff --git a/parameters.cpp b/parameters.cpp index e78365a565..8136aa76e1 100644 --- a/parameters.cpp +++ b/parameters.cpp @@ -71,8 +71,6 @@ Real P::fieldSolverMaxCFL = NAN; Real P::fieldSolverMinCFL = NAN; uint P::fieldSolverSubcycles = 1; -bool P::amrTransShortPencils = false; - uint P::tstep = 0; uint P::tstep_min = 0; uint P::tstep_max = 0; @@ -156,7 +154,8 @@ Realf P::vamrRefineLimit = 1.0; Realf P::vamrCoarsenLimit = 0.5; string P::vamrVelRefCriterion = string(""); -uint P::amrMaxSpatialRefLevel = 0; +bool P::amrTransShortPencils = false; +int P::amrMaxSpatialRefLevel = 0; bool P::adaptRefinement = false; bool P::refineOnRestart = false; bool P::forceRefinement = false; @@ -170,12 +169,14 @@ bool P::useAlpha = false; bool P::useJPerB = false; Real P::JPerBModifier = 0.0; int P::maxFilteringPasses = 0; -uint P::amrBoxHalfWidthX = 1; -uint P::amrBoxHalfWidthY = 1; -uint P::amrBoxHalfWidthZ = 1; -Realf P::amrBoxCenterX = 0.0; -Realf P::amrBoxCenterY = 0.0; -Realf P::amrBoxCenterZ = 0.0; +uint P::amrBoxNumber = 0; +std::vector P::amrBoxHalfWidthX; +std::vector P::amrBoxHalfWidthY; +std::vector P::amrBoxHalfWidthZ; +std::vector P::amrBoxCenterX; +std::vector P::amrBoxCenterY; +std::vector P::amrBoxCenterZ; +std::vector P::amrBoxMaxLevel; vector P::blurPassString; std::vector P::numPasses; //numpasses @@ -446,12 +447,14 @@ bool P::addParameters() { RP::add("AMR.use_alpha","Use alpha as a refinement index", false); RP::add("AMR.use_J_per_B","Use J/B_perp as a refinement index", false); RP::add("AMR.J_per_B_modifier","Factor to add to log2(J / B_perp) in refinement.", 0.0); - RP::add("AMR.box_half_width_x", "Half width of the box that is refined (for testing)", (uint)1); - RP::add("AMR.box_half_width_y", "Half width of the box that is refined (for testing)", (uint)1); - RP::add("AMR.box_half_width_z", "Half width of the box that is refined (for testing)", (uint)1); - RP::add("AMR.box_center_x", "x coordinate of the center of the box that is refined (for testing)", 0.0); - RP::add("AMR.box_center_y", "y coordinate of the center of the box that is refined (for testing)", 0.0); - RP::add("AMR.box_center_z", "z coordinate of the center of the box that is refined (for testing)", 0.0); + RP::add("AMR.number_of_boxes", "How many boxes to be refined, that number of centers and sizes have to then be defined as well.", 0); + RP::addComposing("AMR.box_half_width_x", "Half width in x of the box that is refined"); + RP::addComposing("AMR.box_half_width_y", "Half width in y of the box that is refined"); + RP::addComposing("AMR.box_half_width_z", "Half width in z of the box that is refined"); + RP::addComposing("AMR.box_center_x", "x coordinate of the center of the box that is refined"); + RP::addComposing("AMR.box_center_y", "y coordinate of the center of the box that is refined"); + RP::addComposing("AMR.box_center_z", "z coordinate of the center of the box that is refined"); + RP::addComposing("AMR.box_max_level", "max refinement level of the box that is refined"); RP::add("AMR.transShortPencils", "if true, use one-cell pencils", false); RP::addComposing("AMR.filterpasses", string("AMR filter passes for each individual refinement level")); @@ -690,6 +693,8 @@ void Parameters::getParameters() { RP::get("AMR.use_alpha",P::useAlpha); RP::get("AMR.use_J_per_B",P::useJPerB); RP::get("AMR.J_per_B_modifier",P::JPerBModifier); + RP::get("AMR.number_of_boxes", P::amrBoxNumber); + RP::get("AMR.box_max_level", P::amrBoxMaxLevel); RP::get("AMR.box_half_width_x", P::amrBoxHalfWidthX); RP::get("AMR.box_half_width_y", P::amrBoxHalfWidthY); RP::get("AMR.box_half_width_z", P::amrBoxHalfWidthZ); @@ -699,12 +704,25 @@ void Parameters::getParameters() { RP::get("AMR.transShortPencils", P::amrTransShortPencils); RP::get("AMR.filterpasses", P::blurPassString); + // We need the correct number of parameters for the AMR boxes + if( P::amrBoxNumber != P::amrBoxHalfWidthX.size() + || P::amrBoxNumber != P::amrBoxHalfWidthY.size() + || P::amrBoxNumber != P::amrBoxHalfWidthZ.size() + || P::amrBoxNumber != P::amrBoxCenterX.size() + || P::amrBoxNumber != P::amrBoxCenterY.size() + || P::amrBoxNumber != P::amrBoxCenterZ.size() + || P::amrBoxNumber != P::amrBoxMaxLevel.size() + ) { + cerr << "AMR.number_of_boxes is set to " << P::amrBoxNumber << " so the same number of values is required for AMR.box_half_width_[xyz] and AMR.box_center_[xyz]." << endl; + MPI_Abort(MPI_COMM_WORLD, 1); + } + // If we are in an AMR run we need to set up the filtering scheme. if (P::amrMaxSpatialRefLevel>0){ bool isEmpty = blurPassString.size() == 0; if (!isEmpty){ //sanity check=> user should define a pass for every level - if (blurPassString.size() != P::amrMaxSpatialRefLevel + 1) { + if ((int)blurPassString.size() != P::amrMaxSpatialRefLevel + 1) { cerr << "Filter Passes=" << blurPassString.size() << "\t" << "AMR Levels=" << P::amrMaxSpatialRefLevel + 1 << endl; cerr << "FilterPasses do not match AMR levels. \t" << " in " << __FILE__ << ":" << __LINE__ << endl; MPI_Abort(MPI_COMM_WORLD, 1); @@ -728,7 +746,7 @@ void Parameters::getParameters() { return retval; }; int maxPasses=g_sequence(P::amrMaxSpatialRefLevel-1); - for (uint refLevel=0; refLevel<=P::amrMaxSpatialRefLevel; refLevel++){ + for (int refLevel=0; refLevel<=P::amrMaxSpatialRefLevel; refLevel++){ numPasses.push_back(maxPasses); maxPasses/=2; } diff --git a/parameters.h b/parameters.h index 874098bdfd..4f0658788a 100644 --- a/parameters.h +++ b/parameters.h @@ -185,7 +185,7 @@ struct Parameters { * refined. The value must be larger than vamrCoarsenLimit.*/ static std::string vamrVelRefCriterion; /**< Name of the velocity block refinement criterion function.*/ - static uint amrMaxSpatialRefLevel; + static int amrMaxSpatialRefLevel; static bool adaptRefinement; static bool refineOnRestart; static bool forceRefinement; @@ -199,12 +199,14 @@ struct Parameters { static bool useJPerB; static Real JPerBModifier; static int maxFilteringPasses; - static uint amrBoxHalfWidthX; - static uint amrBoxHalfWidthY; - static uint amrBoxHalfWidthZ; - static Realf amrBoxCenterX; - static Realf amrBoxCenterY; - static Realf amrBoxCenterZ; + static uint amrBoxNumber; + static std::vector amrBoxHalfWidthX; + static std::vector amrBoxHalfWidthY; + static std::vector amrBoxHalfWidthZ; + static std::vector amrBoxCenterX; + static std::vector amrBoxCenterY; + static std::vector amrBoxCenterZ; + static std::vector amrBoxMaxLevel; static bool amrTransShortPencils; /*!< Use short or longpencils in AMR translation.*/ static std::vector blurPassString; static std::vector numPasses; diff --git a/projects/Flowthrough/Flowthrough.cfg b/projects/Flowthrough/Flowthrough.cfg index 146bf111f7..39e39eb864 100644 --- a/projects/Flowthrough/Flowthrough.cfg +++ b/projects/Flowthrough/Flowthrough.cfg @@ -13,9 +13,14 @@ charge = 1 [AMR] max_spatial_level = 2 +number_of_boxes = 1 +box_max_level = 2 box_half_width_x = 2 -box_half_width_z = 2 box_half_width_y = 2 +box_half_width_z = 2 +box_center_x = 0 +box_center_y = 0 +box_center_z = 0 [gridbuilder] x_length = 16 diff --git a/projects/Flowthrough/Flowthrough.cpp b/projects/Flowthrough/Flowthrough.cpp index 1b735d5b9e..58602f3a51 100644 --- a/projects/Flowthrough/Flowthrough.cpp +++ b/projects/Flowthrough/Flowthrough.cpp @@ -266,14 +266,16 @@ namespace projects { std::vector cells = mpiGrid.get_cells(); for (CellID id : cells) { std::array xyz = mpiGrid.get_center(id); - bool inBox = xyz[0] > P::amrBoxCenterX - P::amrBoxHalfWidthX * mpiGrid[id]->parameters[CellParams::DX] && - xyz[0] < P::amrBoxCenterX + P::amrBoxHalfWidthX * mpiGrid[id]->parameters[CellParams::DX] && - xyz[1] > P::amrBoxCenterY - P::amrBoxHalfWidthY * mpiGrid[id]->parameters[CellParams::DY] && - xyz[1] < P::amrBoxCenterY + P::amrBoxHalfWidthY * mpiGrid[id]->parameters[CellParams::DY] && - xyz[2] > P::amrBoxCenterZ - P::amrBoxHalfWidthZ * mpiGrid[id]->parameters[CellParams::DZ] && - xyz[2] < P::amrBoxCenterZ + P::amrBoxHalfWidthZ * mpiGrid[id]->parameters[CellParams::DZ]; - if (inBox) { - refines += mpiGrid.refine_completely(id); + for(uint n = 0; n < P::amrBoxNumber; n++) { + bool inBox = xyz[0] > P::amrBoxCenterX[n] - P::amrBoxHalfWidthX[n] * mpiGrid[id]->parameters[CellParams::DX] && + xyz[0] < P::amrBoxCenterX[n] + P::amrBoxHalfWidthX[n] * mpiGrid[id]->parameters[CellParams::DX] && + xyz[1] > P::amrBoxCenterY[n] - P::amrBoxHalfWidthY[n] * mpiGrid[id]->parameters[CellParams::DY] && + xyz[1] < P::amrBoxCenterY[n] + P::amrBoxHalfWidthY[n] * mpiGrid[id]->parameters[CellParams::DY] && + xyz[2] > P::amrBoxCenterZ[n] - P::amrBoxHalfWidthZ[n] * mpiGrid[id]->parameters[CellParams::DZ] && + xyz[2] < P::amrBoxCenterZ[n] + P::amrBoxHalfWidthZ[n] * mpiGrid[id]->parameters[CellParams::DZ]; + if (inBox) { + refines += mpiGrid.refine_completely(id); + } } } diff --git a/projects/Flowthrough/Flowthrough_amr_test_20190611_YPK.cfg b/projects/Flowthrough/Flowthrough_amr_test_20190611_YPK.cfg index 831a6bbcdc..ae97590f38 100644 --- a/projects/Flowthrough/Flowthrough_amr_test_20190611_YPK.cfg +++ b/projects/Flowthrough/Flowthrough_amr_test_20190611_YPK.cfg @@ -13,9 +13,14 @@ charge = 1 [AMR] max_spatial_level = 1 +number_of_boxes = 1 +box_max_level = 1 +box_center_x = 0 +box_center_y = 0 +box_center_z = 0 box_half_width_x = 1 -box_half_width_z = 1 box_half_width_y = 1 +box_half_width_z = 1 [gridbuilder] x_length = 14 diff --git a/projects/project.cpp b/projects/project.cpp index 287489e8e0..535b354c80 100644 --- a/projects/project.cpp +++ b/projects/project.cpp @@ -524,80 +524,58 @@ namespace projects { cerr << "(Project.cpp) Base class 'refineSpatialCells' in " << __FILE__ << ":" << __LINE__ << " called. Make sure that this is correct." << endl; } - MPI_Comm_rank(MPI_COMM_WORLD,&myRank); - if(myRank == MASTER_RANK) std::cout << "Maximum refinement level is " << mpiGrid.mapping.get_maximum_refinement_level() << std::endl; std::vector refineSuccess; - for (uint i = 0; i < 2 * P::amrBoxHalfWidthX; ++i) { - for (uint j = 0; j < 2 * P::amrBoxHalfWidthY; ++j) { - for (uint k = 0; k < 2 * P::amrBoxHalfWidthZ; ++k) { - - std::array xyz; - xyz[0] = P::amrBoxCenterX + (0.5 + i - P::amrBoxHalfWidthX) * P::dx_ini; - xyz[1] = P::amrBoxCenterY + (0.5 + j - P::amrBoxHalfWidthY) * P::dy_ini; - xyz[2] = P::amrBoxCenterZ + (0.5 + k - P::amrBoxHalfWidthZ) * P::dz_ini; - - if (mpiGrid.refine_completely_at(xyz)) { - #ifndef NDEBUG - CellID myCell = mpiGrid.get_existing_cell(xyz); - std::cout << "Rank " << myRank << " is refining cell " << myCell << std::endl; - #endif + for (int level = 0; level < mpiGrid.mapping.get_maximum_refinement_level(); level++) { + int refineCount = 0; + for (int n = 0; n < P::amrBoxNumber; n++) { + if (level < P::amrBoxMaxLevel[n]) { + for (int i = 0; i < pow(2, level+1) * P::amrBoxHalfWidthX[n]; ++i) { + for (int j = 0; j < pow(2, level+1) * P::amrBoxHalfWidthY[n]; ++j) { + for (int k = 0; k < pow(2, level+1) * P::amrBoxHalfWidthZ[n]; ++k) { + + std::array xyz; + xyz[0] = P::amrBoxCenterX[n] + (0.5 + i - pow(2, level)*P::amrBoxHalfWidthX[n]) * P::dx_ini / pow(2, level); + xyz[1] = P::amrBoxCenterY[n] + (0.5 + j - pow(2, level)*P::amrBoxHalfWidthY[n]) * P::dy_ini / pow(2, level); + xyz[2] = P::amrBoxCenterZ[n] + (0.5 + k - pow(2, level)*P::amrBoxHalfWidthZ[n]) * P::dz_ini / pow(2, level); + + if (mpiGrid.refine_completely_at(xyz)) { + refineCount++; + #ifndef NDEBUG + CellID myCell = mpiGrid.get_existing_cell(xyz); + std::cout << "Rank " << myRank << " is refining cell " << myCell << std::endl; + #endif + } // if + } // box z + } // box y + } // box x + } // if (P::amrBoxMaxLevel <= level) + } // box number + int totalRefineCount; + MPI_Allreduce(&refineCount, &totalRefineCount, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + if(totalRefineCount > 0) { + std::vector refinedCells = mpiGrid.stop_refining(true); + + #ifndef NDEBUG + if(refinedCells.size() > 0) { + std::cerr << "Refined cells produced by rank " << myRank << " for level " << level << " are: "; + for (auto cellid : refinedCells) { + std::cout << cellid << " "; } + std::cout << endl; } + #endif + + mpiGrid.balance_load(); } - } - std::vector refinedCells = mpiGrid.stop_refining(true); - if(myRank == MASTER_RANK) std::cout << "Finished first level of refinement" << endl; - #ifndef NDEBUG - if(refinedCells.size() > 0) { - std::cout << "Refined cells produced by rank " << myRank << " are: "; - for (auto cellid : refinedCells) { - std::cout << cellid << " "; + if(myRank == MASTER_RANK) { + std::cout << "Finished level of refinement " << level+1 << endl; } - std::cout << endl; - } - #endif - - mpiGrid.balance_load(); - - if(mpiGrid.get_maximum_refinement_level() > 1) { - for (uint i = 0; i < 2 * P::amrBoxHalfWidthX; ++i) { - for (uint j = 0; j < 2 * P::amrBoxHalfWidthY; ++j) { - for (uint k = 0; k < 2 * P::amrBoxHalfWidthZ; ++k) { - - std::array xyz; - xyz[0] = P::amrBoxCenterX + 0.5 * (0.5 + i - P::amrBoxHalfWidthX) * P::dx_ini; - xyz[1] = P::amrBoxCenterY + 0.5 * (0.5 + j - P::amrBoxHalfWidthY) * P::dy_ini; - xyz[2] = P::amrBoxCenterZ + 0.5 * (0.5 + k - P::amrBoxHalfWidthZ) * P::dz_ini; - - if (mpiGrid.refine_completely_at(xyz)) { - #ifndef NDEBUG - CellID myCell = mpiGrid.get_existing_cell(xyz); - std::cout << "Rank " << myRank << " is refining cell " << myCell << std::endl; - #endif - } - } - } - } - - std::vector refinedCells = mpiGrid.stop_refining(true); - if(myRank == MASTER_RANK) std::cout << "Finished second level of refinement" << endl; - #ifndef NDEBUG - if(refinedCells.size() > 0) { - std::cout << "Refined cells produced by rank " << myRank << " are: "; - for (auto cellid : refinedCells) { - std::cout << cellid << " "; - } - std::cout << endl; - } - #endif - mpiGrid.balance_load(); - } - - return true; + } // refinement levels + return true; } int Project::adaptRefinement( dccrg::Dccrg& mpiGrid ) const { diff --git a/projects/testAmr/testAmr.cfg b/projects/testAmr/testAmr.cfg index 8b52162d9a..6a83ee3c89 100644 --- a/projects/testAmr/testAmr.cfg +++ b/projects/testAmr/testAmr.cfg @@ -23,6 +23,8 @@ system_write_distribution_zline_stride = 0 [AMR] max_spatial_level = 2 +number_of_boxes = 1 +box_max_level = 1 box_half_width_x = 1 box_half_width_y = 1 box_half_width_z = 1 @@ -100,4 +102,4 @@ rhoPertAbsAmp = 0.0 [loadBalance] #algorithm = RCB -algorithm = RANDOM \ No newline at end of file +algorithm = RANDOM diff --git a/projects/testAmr/testAmr.cpp b/projects/testAmr/testAmr.cpp index e873b0450e..9508b95f85 100644 --- a/projects/testAmr/testAmr.cpp +++ b/projects/testAmr/testAmr.cpp @@ -62,7 +62,6 @@ namespace projects { RP::add("testAmr.lambda", "B cosine perturbation wavelength (m)", 1.0); RP::add("testAmr.nVelocitySamples", "Number of sampling points per velocity dimension", 2); RP::add("testAmr.densityModel","Which spatial density model is used?",string("uniform")); - RP::add("testAmr.maxSpatialRefinementLevel", "Maximum level for spatial refinement", 1.0); // Per-population parameters for(uint i=0; i< getObjectWrapper().particleSpecies.size(); i++) { @@ -93,7 +92,6 @@ namespace projects { RP::get("testAmr.dBy", this->dBy); RP::get("testAmr.dBz", this->dBz); RP::get("testAmr.lambda", this->lambda); - RP::get("testAmr.maxSpatialRefinementLevel", this->maxSpatialRefinementLevel); RP::get("testAmr.nVelocitySamples", this->nVelocitySamples); // Per-population parameters @@ -273,84 +271,5 @@ namespace projects { } return centerPoints; } - - bool testAmr::refineSpatialCells( dccrg::Dccrg& mpiGrid ) const { - - int myRank; - MPI_Comm_rank(MPI_COMM_WORLD,&myRank); - - if(myRank == MASTER_RANK) std::cout << "Maximum refinement level is " << mpiGrid.mapping.get_maximum_refinement_level() << std::endl; - - std::vector refineSuccess; - - for (uint i = 0; i < 2 * P::amrBoxHalfWidthX; ++i) { - for (uint j = 0; j < 2 * P::amrBoxHalfWidthY; ++j) { - for (uint k = 0; k < 2 * P::amrBoxHalfWidthZ; ++k) { - - std::array xyz; - xyz[0] = P::amrBoxCenterX + (0.5 + i - P::amrBoxHalfWidthX) * P::dx_ini; - xyz[1] = P::amrBoxCenterY + (0.5 + j - P::amrBoxHalfWidthY) * P::dy_ini; - xyz[2] = P::amrBoxCenterZ + (0.5 + k - P::amrBoxHalfWidthZ) * P::dz_ini; - - if (mpiGrid.refine_completely_at(xyz)) { -#ifndef NDEBUG - CellID myCell = mpiGrid.get_existing_cell(xyz); - std::cout << "Rank " << myRank << " is refining cell " << myCell << std::endl; -#endif - } - } - } - } - std::vector refinedCells = mpiGrid.stop_refining(true); - if(myRank == MASTER_RANK) std::cout << "Finished first level of refinement" << endl; -#ifndef NDEBUG - if(refinedCells.size() > 0) { - std::cout << "Refined cells produced by rank " << myRank << " are: "; - for (auto cellid : refinedCells) { - std::cout << cellid << " "; - } - std::cout << endl; - } -#endif - - mpiGrid.balance_load(); - - if(mpiGrid.get_maximum_refinement_level() > 1) { - - for (uint i = 0; i < 2 * P::amrBoxHalfWidthX; ++i) { - for (uint j = 0; j < 2 * P::amrBoxHalfWidthY; ++j) { - for (uint k = 0; k < 2 * P::amrBoxHalfWidthZ; ++k) { - - std::array xyz; - xyz[0] = P::amrBoxCenterX + 0.5 * (0.5 + i - P::amrBoxHalfWidthX) * P::dx_ini; - xyz[1] = P::amrBoxCenterY + 0.5 * (0.5 + j - P::amrBoxHalfWidthY) * P::dy_ini; - xyz[2] = P::amrBoxCenterZ + 0.5 * (0.5 + k - P::amrBoxHalfWidthZ) * P::dz_ini; - - if (mpiGrid.refine_completely_at(xyz)) { -#ifndef NDEBUG - CellID myCell = mpiGrid.get_existing_cell(xyz); - std::cout << "Rank " << myRank << " is refining cell " << myCell << std::endl; -#endif - } - } - } - } - - std::vector refinedCells = mpiGrid.stop_refining(true); - if(myRank == MASTER_RANK) std::cout << "Finished second level of refinement" << endl; -#ifndef NDEBUG - if(refinedCells.size() > 0) { - std::cout << "Refined cells produced by rank " << myRank << " are: "; - for (auto cellid : refinedCells) { - std::cout << cellid << " "; - } - std::cout << endl; - } -#endif - mpiGrid.balance_load(); - } - - return true; - } }// namespace projects diff --git a/projects/testAmr/testAmr.h b/projects/testAmr/testAmr.h index 51e295ef7d..3fbbabd96d 100644 --- a/projects/testAmr/testAmr.h +++ b/projects/testAmr/testAmr.h @@ -84,7 +84,6 @@ namespace projects { creal z, const uint popID ) const; - bool refineSpatialCells( dccrg::Dccrg& mpiGrid ) const; static Real rhoRnd; //static as it has to be threadprivate #pragma omp threadprivate(rhoRnd) @@ -98,7 +97,6 @@ namespace projects { Real magYPertAbsAmp; Real magZPertAbsAmp; Real lambda; - int maxSpatialRefinementLevel; uint nVelocitySamples; std::vector speciesParams; diff --git a/projects/test_fp/test_fp.cpp b/projects/test_fp/test_fp.cpp index 73efe181b9..5476a60031 100644 --- a/projects/test_fp/test_fp.cpp +++ b/projects/test_fp/test_fp.cpp @@ -274,57 +274,4 @@ namespace projects { return this->getV0(x,y,z,dx,dy,dz,popID); } - bool test_fp::refineSpatialCells( dccrg::Dccrg& mpiGrid ) const { - - int myRank; - MPI_Comm_rank(MPI_COMM_WORLD,&myRank); - - // mpiGrid.set_maximum_refinement_level(std::min(this->maxSpatialRefinementLevel, mpiGrid.mapping.get_maximum_refinement_level())); - - // cout << "I am at line " << __LINE__ << " of " << __FILE__ << endl; - if(myRank == MASTER_RANK) std::cout << "Maximum refinement level is " << mpiGrid.mapping.get_maximum_refinement_level() << std::endl; - - - for (double x = P::amrBoxCenterX - P::amrBoxHalfWidthX * P::dx_ini; x <= P::amrBoxCenterX + P::amrBoxHalfWidthX * P::dx_ini; x += 0.99 * P::dx_ini) { - for (double y = P::amrBoxCenterY - P::amrBoxHalfWidthY * P::dy_ini; y <= P::amrBoxCenterY + P::amrBoxHalfWidthY * P::dy_ini; y += 0.99 * P::dy_ini) { - for (double z = P::amrBoxCenterZ - P::amrBoxHalfWidthZ * P::dz_ini; z <= P::amrBoxCenterZ + P::amrBoxHalfWidthZ * P::dz_ini; z += 0.99 * P::dz_ini) { - - std::array xyz; - xyz[0] = x; - xyz[1] = y; - xyz[2] = z; - CellID myCell = mpiGrid.get_existing_cell(xyz); - if (mpiGrid.refine_completely_at(xyz)) { - std::cout << "Rank " << myRank << " is refining cell " << myCell << std::endl; - } - } - } - } - - std::vector refinedCells = mpiGrid.stop_refining(true); - if(myRank == MASTER_RANK) std::cout << "Finished first level of refinement" << endl; - if(refinedCells.size() > 0) { - std::cout << "Refined cells produced by rank " << myRank << " are: "; - for (auto cellid : refinedCells) { - std::cout << cellid << " "; - } - std::cout << endl; - } - - mpiGrid.balance_load(); - -// const vector& cells = getLocalCells(); -// if(cells.empty()) { -// std::cout << "Rank " << myRank << " has no cells!" << std::endl; -// } else { -// std::cout << "Cells on rank " << myRank << ": "; -// for (auto c : cells) { -// std::cout << c << " "; -// } -// std::cout << std::endl; -// } - - return true; - } - }// namespace projects diff --git a/projects/test_fp/test_fp.h b/projects/test_fp/test_fp.h index 22adb53f68..af1e4bbb19 100644 --- a/projects/test_fp/test_fp.h +++ b/projects/test_fp/test_fp.h @@ -46,7 +46,6 @@ namespace projects { protected: Real sign(creal value) const; Real getDistribValue(creal& vx, creal& vy, creal& vz); - bool refineSpatialCells( dccrg::Dccrg& mpiGrid ) const; virtual void calcCellParameters(spatial_cell::SpatialCell* cell,creal& t); virtual Real calcPhaseSpaceDensity( creal& x, creal& y, creal& z, diff --git a/spatial_cell.cpp b/spatial_cell.cpp index 8efa634d1e..d1ddd47dfb 100644 --- a/spatial_cell.cpp +++ b/spatial_cell.cpp @@ -605,8 +605,13 @@ namespace spatial_cell { } // send velocity block list - displacements.push_back((uint8_t*) &(populations[activePopID].vmesh.getGrid()[0]) - (uint8_t*) this); - block_lengths.push_back(sizeof(vmesh::GlobalID) * populations[activePopID].vmesh.size()); + if(populations[activePopID].vmesh.size() > 0) { + displacements.push_back((uint8_t*) &(populations[activePopID].vmesh.getGrid()[0]) - (uint8_t*) this); + block_lengths.push_back(sizeof(vmesh::GlobalID) * populations[activePopID].vmesh.size()); + } else { + displacements.push_back(0); + block_lengths.push_back(0); + } } if ((SpatialCell::mpi_transfer_type & Transfer::VEL_BLOCK_WITH_CONTENT_STAGE1) !=0) { @@ -621,8 +626,13 @@ namespace spatial_cell { } //velocity_block_with_content_list_size should first be updated, before this can be done (STAGE1) - displacements.push_back((uint8_t*) &(this->velocity_block_with_content_list[0]) - (uint8_t*) this); - block_lengths.push_back(sizeof(vmesh::GlobalID)*this->velocity_block_with_content_list_size); + if(velocity_block_with_content_list_size > 0) { + displacements.push_back((uint8_t*) &(this->velocity_block_with_content_list[0]) - (uint8_t*) this); + block_lengths.push_back(sizeof(vmesh::GlobalID)*this->velocity_block_with_content_list_size); + } else { + displacements.push_back(0); + block_lengths.push_back(0); + } } if ((SpatialCell::mpi_transfer_type & Transfer::VEL_BLOCK_DATA) !=0) { @@ -642,6 +652,13 @@ namespace spatial_cell { if ( P::amrMaxSpatialRefLevel == 0 || receiving || ranks.find(receiver_rank) != ranks.end()) { for ( int i = 0; i < MAX_NEIGHBORS_PER_DIM; ++i) { + if(!receiving) { + for(unsigned int j=0; j < this->neighbor_number_of_blocks[i]*VELOCITY_BLOCK_LENGTH; j++) { + if(isnan(this->neighbor_block_data[i][j])) { + std::cerr << "WARNING: Sending NaN data to neighbour (cellID " << cellID << ", neighbourID " << i << ", neigbour block data " << j << std::endl; + } + } + } displacements.push_back((uint8_t*) this->neighbor_block_data[i] - (uint8_t*) this); block_lengths.push_back(sizeof(Realf) * VELOCITY_BLOCK_LENGTH * this->neighbor_number_of_blocks[i]); } diff --git a/testpackage/tests/Flowthrough_amr/Flowthrough_amr.cfg b/testpackage/tests/Flowthrough_amr/Flowthrough_amr.cfg index fa7cf0f4eb..5d5169bd96 100644 --- a/testpackage/tests/Flowthrough_amr/Flowthrough_amr.cfg +++ b/testpackage/tests/Flowthrough_amr/Flowthrough_amr.cfg @@ -13,9 +13,14 @@ charge = 1 [AMR] max_spatial_level = 2 +number_of_boxes = 1 +box_max_level = 2 box_half_width_x = 1 -box_half_width_z = 1 box_half_width_y = 1 +box_half_width_z = 1 +box_center_x = 0 +box_center_y = 0 +box_center_z = 0 [gridbuilder] x_length = 16 diff --git a/testpackage/tests/Flowthrough_trans_periodic/Flowthrough_trans_periodic.cfg b/testpackage/tests/Flowthrough_trans_periodic/Flowthrough_trans_periodic.cfg index fea58d0daa..66685b549c 100644 --- a/testpackage/tests/Flowthrough_trans_periodic/Flowthrough_trans_periodic.cfg +++ b/testpackage/tests/Flowthrough_trans_periodic/Flowthrough_trans_periodic.cfg @@ -9,7 +9,7 @@ ParticlePopulations = proton [io] write_initial_state = 0 -system_write_t_interval = 649.0 +system_write_t_interval = 39.44 system_write_file_name = bulk system_write_distribution_stride = 1 system_write_distribution_xline_stride = 0 @@ -30,15 +30,15 @@ output = populations_vg_nonmaxwellianity diagnostic = populations_vg_blocks [gridbuilder] -x_length = 20 -y_length = 20 +x_length = 8 +y_length = 1 z_length = 1 -x_min = -1.3e8 -x_max = 1.3e8 -y_min = -1.3e8 -y_max = 1.3e8 -z_min = -6.5e6 -z_max = 6.5e6 +x_min = -10.4e7 +x_max = 10.4e7 +y_min = -1.3e7 +y_max = 1.3e7 +z_min = -1.3e7 +z_max = 1.3e7 t_max = 650 dt = 2.0 @@ -54,9 +54,9 @@ vy_min = -600000.0 vy_max = +600000.0 vz_min = -600000.0 vz_max = +600000.0 -vx_length = 15 -vy_length = 15 -vz_length = 15 +vx_length = 5 +vy_length = 5 +vz_length = 5 [proton_sparse] minValue = 1.0e-15 @@ -78,8 +78,8 @@ densityModel = SheetMaxwellian T = 100000.0 rho = 1000000.0 VX0 = 4e5 -VY0 = 4e5 -VZ0 = 4e5 +VY0 = 0 +VZ0 = 0 nSpaceSamples = 2 nVelocitySamples = 2 diff --git a/testpackage/tests/Magnetosphere_small/Magnetosphere_small.cfg b/testpackage/tests/Magnetosphere_small/Magnetosphere_small.cfg index c820e044b6..5243bc7abc 100644 --- a/testpackage/tests/Magnetosphere_small/Magnetosphere_small.cfg +++ b/testpackage/tests/Magnetosphere_small/Magnetosphere_small.cfg @@ -11,7 +11,7 @@ charge = 1 diagnostic_write_interval = 1 write_initial_state = 0 -system_write_t_interval = 10 +system_write_t_interval = .027 system_write_file_name = bulk system_write_distribution_stride = 0 system_write_distribution_xline_stride = 10 diff --git a/testpackage/tests/transtest_amr/transtest_amr.cfg b/testpackage/tests/transtest_amr/transtest_amr.cfg index bb065566ac..eead3b0ed0 100644 --- a/testpackage/tests/transtest_amr/transtest_amr.cfg +++ b/testpackage/tests/transtest_amr/transtest_amr.cfg @@ -23,6 +23,8 @@ system_write_distribution_zline_stride = 0 [AMR] max_spatial_level = 2 +number of boxes = 1 +box_max_level = 2 box_half_width_x = 1 box_half_width_y = 1 box_half_width_z = 1 diff --git a/vlasiator.cpp b/vlasiator.cpp index fae1a394ae..7517a49c09 100644 --- a/vlasiator.cpp +++ b/vlasiator.cpp @@ -483,7 +483,7 @@ int main(int argn,char* args[]) { sysBoundaryContainer, *project ); - + const std::vector& cells = getLocalCells(); initGridsTimer.stop(); @@ -747,7 +747,7 @@ int main(int argn,char* args[]) { //report filtering if we are in an AMR run if (P::amrMaxSpatialRefLevel>0){ logFile<<"Filtering Report: "< Passes "<::indices_t indices_unsigned = mpiGrid.mapping.get_indices(cellID); - int64_t indices[3]; - dccrg::Grid_Length::type length = mpiGrid.mapping.length.get(); - - //compute raw new indices - indices[0] = spatial_di + indices_unsigned[0]; - indices[1] = spatial_dj + indices_unsigned[1]; - indices[2] = spatial_dk + indices_unsigned[2]; - - //take periodicity into account - for(uint i = 0; i<3; i++) { - if(mpiGrid.topology.is_periodic(i)) { - while(indices[i] < 0 ) - indices[i] += length[i]; - while(indices[i] >= static_cast(length[i]) ) - indices[i] -= length[i]; + dccrg::Types<3>::indices_t indices = mpiGrid.mapping.get_indices(cellID); + + std::set nbrIDs = mpiGrid.find_cells_at_offset(indices, cellID, 0, {spatial_di, spatial_dj, spatial_dk}); + + if(nbrIDs.size() != 1) { + std::cerr << "Error: Cell " << cellID << " has more than one neighbour (namely " << nbrIDs.size() << ":" << std::endl; + std::cerr << "["; + for(auto n : nbrIDs) { + std::cerr << n << ", "; } + std::cerr << "]" << std::endl; + + abort(); } - //return INVALID_CELLID for cells outside system (non-periodic) - for(uint i = 0; i<3; i++) { - if(indices[i]< 0) - return INVALID_CELLID; - if(indices[i]>=static_cast(length[i])) - return INVALID_CELLID; - } - //store nbr indices into the correct datatype - for(uint i = 0; i<3; i++) { - indices_unsigned[i] = indices[i]; - } + //get nbrID - CellID nbrID = mpiGrid.mapping.get_cell_from_indices(indices_unsigned,0); + CellID nbrID = *nbrIDs.begin(); if (nbrID == dccrg::error_cell ) { - std::cerr << __FILE__ << ":" << __LINE__ - << " No neighbor for cell?" << cellID - << " at offsets " << spatial_di << ", " << spatial_dj << ", " << spatial_dk - << std::endl; - abort(); + //std::cerr << __FILE__ << ":" << __LINE__ + // << " No neighbor for cell " << cellID << " (indices [" << indices[0] << ", " << indices[1] << ", " << indices[2] << "]" + // << " at offsets [" << spatial_di << ", " << spatial_dj << ", " << spatial_dk + // << "]" + // << std::endl; + //abort(); + return INVALID_CELLID; } // not existing cell or do not compute - if( mpiGrid[nbrID]->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) + if( mpiGrid[nbrID]->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { + //std::cerr << "Cell " << cellID << " has no neighbour in offset [" << spatial_di << ", " << spatial_dj << ", " << spatial_dk << "] " + // "because the cell would be DO_NOT_COMPUTE" << std::endl; return INVALID_CELLID; + } //cell on boundary, but not first layer and we want to include //first layer (e.g. when we compute source cells) if( include_first_boundary_layer && mpiGrid[nbrID]->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && mpiGrid[nbrID]->sysBoundaryLayer != 1 ) { + //std::cerr << "Cell " << cellID << " has no neighbour in offset [" << spatial_di << ", " << spatial_dj << ", " << spatial_dk << "] " + // "because the cell would be boundary layer 1." << std::endl; return INVALID_CELLID; } @@ -132,9 +125,12 @@ CellID get_spatial_neighbor(const dccrg::DccrgsysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY){ + //std::cerr << "Cell " << cellID << " has no neighbour in offset [" << spatial_di << ", " << spatial_dj << ", " << spatial_dk << "] " + // "because the cell would be on the boundary." << std::endl; return INVALID_CELLID; } + //std::cerr << "Cell " << cellID << " has neighbour " << nbrID << "in offset [" << spatial_di << ", " << spatial_dj << ", " << spatial_dk << "] " << std::endl; return nbrID; //no AMR } @@ -682,23 +678,44 @@ void update_remote_mapping_contribution( CellID p_ngbr = INVALID_CELLID; CellID m_ngbr = INVALID_CELLID; + //if(local_cells[c] == 4 || local_cells[c] == 5) { + // fprintf(stderr, "Cell %li has neighbours in SHIFT_P_X_NEIGHBORHOOD_ID: [\n", local_cells[c]); + // + // for(auto& nbr : *mpiGrid.get_neighbors_of(local_cells[c], SHIFT_P_X_NEIGHBORHOOD_ID)) { + // fprintf(stderr, "\t%li (%i %i %i %i)\n", nbr.first, nbr.second[0], nbr.second[1], nbr.second[2], nbr.second[3]); + // } + // fprintf(stderr, "]\n"); + // fprintf(stderr, "Cell %li has neighbours in SHIFT_M_X_NEIGHBORHOOD_ID: [\n", local_cells[c]); + // for(auto& nbr : *mpiGrid.get_neighbors_of(local_cells[c], SHIFT_M_X_NEIGHBORHOOD_ID)) { + // fprintf(stderr, "\t%li (%i %i %i %i)\n", nbr.first, nbr.second[0], nbr.second[1], nbr.second[2], nbr.second[3]); + // } + // fprintf(stderr, "]\n"); + + //} + for (const auto& [neighbor, dir] : mpiGrid.get_face_neighbors_of(local_cells[c])) { if(dir == ((int)dimension + 1) * direction) { p_ngbr = neighbor; } - if(dir == -1 * ((int)dimension + 1) * direction) { m_ngbr = neighbor; } + } - + + //MPI_Barrier(MPI_COMM_WORLD); + //internal cell, not much to do if (mpiGrid.is_local(p_ngbr) && mpiGrid.is_local(m_ngbr)) continue; SpatialCell *pcell = NULL; - if (p_ngbr != INVALID_CELLID) pcell = mpiGrid[p_ngbr]; + if (p_ngbr != INVALID_CELLID) { + pcell = mpiGrid[p_ngbr]; + } SpatialCell *mcell = NULL; - if (m_ngbr != INVALID_CELLID) mcell = mpiGrid[m_ngbr]; + if (m_ngbr != INVALID_CELLID) { + mcell = mpiGrid[m_ngbr]; + } if (p_ngbr != INVALID_CELLID && pcell->sysBoundaryFlag == sysboundarytype::NOT_SYSBOUNDARY) if (!mpiGrid.is_local(p_ngbr) && do_translate_cell(ccell)) { //if (p_ngbr != INVALID_CELLID && !mpiGrid.is_local(p_ngbr) && do_translate_cell(ccell)) { @@ -707,12 +724,22 @@ void update_remote_mapping_contribution( //2) is remote cell, 3) if the source cell in center was //translated ccell->neighbor_block_data[0] = pcell->get_data(popID); + + for(unsigned int cell = 0; cellget_number_of_velocity_blocks(popID); ++cell) { + if(isnan( pcell->get_data(popID)[cell] ) || isinf( pcell->get_data(popID)[cell])) { + fprintf(stderr,"NaN sent at cell %li, vel cell %i",receive_cells[c], cell); + abort(); + } + } + + ccell->neighbor_number_of_blocks[0] = pcell->get_number_of_velocity_blocks(popID); send_cells.push_back(p_ngbr); } if (m_ngbr != INVALID_CELLID && !mpiGrid.is_local(m_ngbr) && ccell->sysBoundaryFlag == sysboundarytype::NOT_SYSBOUNDARY) { + //Receive data that mcell mapped to ccell to this local cell //data array, if 1) m is a valid source cell, 2) center cell is to be updated (normal cell) 3) m is remote //we will here allocate a receive buffer, since we need to aggregate values @@ -721,9 +748,11 @@ void update_remote_mapping_contribution( receive_cells.push_back(local_cells[c]); receiveBuffers.push_back(mcell->neighbor_block_data[0]); - } + } } + //std::cerr << "send_cells.size = " << send_cells.size() << ", recieve_cells.size = " << receive_cells.size() << std::endl; + // Do communication SpatialCell::setCommunicatedSpecies(popID); SpatialCell::set_mpi_transfer_type(Transfer::NEIGHBOR_VEL_BLOCK_DATA); @@ -742,16 +771,30 @@ void update_remote_mapping_contribution( break; } -#pragma omp parallel { //reduce data: sum received data in the data array to // the target grid in the temporary block container + //std::cerr << "Recieve_cells are [ "; + //for (size_t c=0; c < receive_cells.size(); ++c) { + // std::cerr << receive_cells[c] << " "; + //} + //std::cerr << "\n"; for (size_t c=0; c < receive_cells.size(); ++c) { SpatialCell* spatial_cell = mpiGrid[receive_cells[c]]; Realf *blockData = spatial_cell->get_data(popID); + //fprintf(stderr, "Cell %i's recieveBuffer is at 0x%08lx, recieveBuffers[c][cell] is at 0x%08lx\n", + // receive_cells[c], + // mpiGrid[(( receive_cells[c]+((direction>0)?-1:1) - 1 ) % 8 ) + 1 ]->neighbor_block_data[0], + // receiveBuffers[c]); + #pragma omp for for(unsigned int cell = 0; cellget_number_of_velocity_blocks(popID); ++cell) { + if(isnan(receiveBuffers[c][cell]) || isinf(receiveBuffers[c][cell])) { + fprintf(stderr, "NaN received at cell %li, vel_cell %i (%i blocks)\n", receive_cells[c], cell, spatial_cell->get_number_of_velocity_blocks(popID)); + //abort(); + //break; + } blockData[cell] += receiveBuffers[c][cell]; } } diff --git a/vlasovsolver/vlasovmover.cpp b/vlasovsolver/vlasovmover.cpp index beade08e51..987c73ad7e 100644 --- a/vlasovsolver/vlasovmover.cpp +++ b/vlasovsolver/vlasovmover.cpp @@ -90,45 +90,45 @@ void calculateSpatialTranslation( btzTimer.stop(); // ------------- SLICE - map dist function in Z --------------- // - if(P::zcells_ini > 1){ + //if(P::zcells_ini > 1){ phiprof::Timer transTimer {"transfer-stencil-data-z", {"MPI"}}; //updateRemoteVelocityBlockLists(mpiGrid,popID,VLASOV_SOLVER_Z_NEIGHBORHOOD_ID); - SpatialCell::set_mpi_transfer_direction(2); - SpatialCell::set_mpi_transfer_type(Transfer::VEL_BLOCK_DATA,false,AMRtranslationActive); - mpiGrid.update_copies_of_remote_neighbors(VLASOV_SOLVER_Z_NEIGHBORHOOD_ID); + //SpatialCell::set_mpi_transfer_direction(2); + //SpatialCell::set_mpi_transfer_type(Transfer::VEL_BLOCK_DATA,false,AMRtranslationActive); + //mpiGrid.update_copies_of_remote_neighbors(VLASOV_SOLVER_Z_NEIGHBORHOOD_ID); transTimer.stop(); - // bt=phiprof::initializeTimer("barrier-trans-pre-trans_map_1d-z","Barriers","MPI"); - // phiprof::start(bt); - // MPI_Barrier(MPI_COMM_WORLD); - // phiprof::stop(bt); + // // bt=phiprof::initializeTimer("barrier-trans-pre-trans_map_1d-z","Barriers","MPI"); + // // phiprof::start(bt); + // // MPI_Barrier(MPI_COMM_WORLD); + // // phiprof::stop(bt); t1 = MPI_Wtime(); phiprof::Timer computeTimer {"compute-mapping-z"}; - if(P::amrMaxSpatialRefLevel == 0) { - trans_map_1d(mpiGrid,local_propagated_cells, remoteTargetCellsz, 2, dt,popID); // map along z// - } else { - trans_map_1d_amr(mpiGrid,local_propagated_cells, remoteTargetCellsz, nPencils, 2, dt,popID); // map along z// - } + //if(P::amrMaxSpatialRefLevel == 0) { + // trans_map_1d(mpiGrid,local_propagated_cells, remoteTargetCellsz, 2, dt,popID); // map along z// + //} else { + // trans_map_1d_amr(mpiGrid,local_propagated_cells, remoteTargetCellsz, nPencils, 2, dt,popID); // map along z// + //} computeTimer.stop(); time += MPI_Wtime() - t1; phiprof::Timer btTimer {"barrier-trans-pre-update_remote-z", {"Barriers","MPI"}}; - MPI_Barrier(MPI_COMM_WORLD); + //MPI_Barrier(MPI_COMM_WORLD); btTimer.stop(); phiprof::Timer updateRemoteTimer {"update_remote-z", {"MPI"}}; - if(P::amrMaxSpatialRefLevel == 0) { - update_remote_mapping_contribution(mpiGrid, 2,+1,popID); - update_remote_mapping_contribution(mpiGrid, 2,-1,popID); - } else { - update_remote_mapping_contribution_amr(mpiGrid, 2,+1,popID); - update_remote_mapping_contribution_amr(mpiGrid, 2,-1,popID); - } + //if(P::amrMaxSpatialRefLevel == 0) { + // update_remote_mapping_contribution(mpiGrid, 2,+1,popID); + // update_remote_mapping_contribution(mpiGrid, 2,-1,popID); + //} else { + // update_remote_mapping_contribution_amr(mpiGrid, 2,+1,popID); + // update_remote_mapping_contribution_amr(mpiGrid, 2,-1,popID); + //} updateRemoteTimer.stop(); - } + //} phiprof::Timer btxTimer {"barrier-trans-pre-x", {"Barriers","MPI"}}; MPI_Barrier(MPI_COMM_WORLD); @@ -183,9 +183,9 @@ void calculateSpatialTranslation( phiprof::Timer transTimer {"transfer-stencil-data-y", {"MPI"}}; //updateRemoteVelocityBlockLists(mpiGrid,popID,VLASOV_SOLVER_Y_NEIGHBORHOOD_ID); - SpatialCell::set_mpi_transfer_direction(1); - SpatialCell::set_mpi_transfer_type(Transfer::VEL_BLOCK_DATA,false,AMRtranslationActive); - mpiGrid.update_copies_of_remote_neighbors(VLASOV_SOLVER_Y_NEIGHBORHOOD_ID); + //SpatialCell::set_mpi_transfer_direction(1); + //SpatialCell::set_mpi_transfer_type(Transfer::VEL_BLOCK_DATA,false,AMRtranslationActive); + //mpiGrid.update_copies_of_remote_neighbors(VLASOV_SOLVER_Y_NEIGHBORHOOD_ID); transTimer.stop(); // bt=phiprof::initializeTimer("barrier-trans-pre-trans_map_1d-y","Barriers","MPI"); @@ -195,26 +195,26 @@ void calculateSpatialTranslation( t1 = MPI_Wtime(); phiprof::Timer computeTimer {"compute-mapping-y"}; - if(P::amrMaxSpatialRefLevel == 0) { - trans_map_1d(mpiGrid,local_propagated_cells, remoteTargetCellsy, 1,dt,popID); // map along y// - } else { - trans_map_1d_amr(mpiGrid,local_propagated_cells, remoteTargetCellsy, nPencils, 1,dt,popID); // map along y// - } + //if(P::amrMaxSpatialRefLevel == 0) { + // trans_map_1d(mpiGrid,local_propagated_cells, remoteTargetCellsy, 1,dt,popID); // map along y// + //} else { + // trans_map_1d_amr(mpiGrid,local_propagated_cells, remoteTargetCellsy, nPencils, 1,dt,popID); // map along y// + //} computeTimer.stop(); time += MPI_Wtime() - t1; phiprof::Timer btTimer {"barrier-trans-pre-update_remote-y", {"Barriers","MPI"}}; - MPI_Barrier(MPI_COMM_WORLD); + //MPI_Barrier(MPI_COMM_WORLD); btTimer.stop(); phiprof::Timer updateRemoteTimer {"update_remote-y", {"MPI"}}; - if(P::amrMaxSpatialRefLevel == 0) { - update_remote_mapping_contribution(mpiGrid, 1,+1,popID); - update_remote_mapping_contribution(mpiGrid, 1,-1,popID); - } else { - update_remote_mapping_contribution_amr(mpiGrid, 1,+1,popID); - update_remote_mapping_contribution_amr(mpiGrid, 1,-1,popID); - } + //if(P::amrMaxSpatialRefLevel == 0) { + // update_remote_mapping_contribution(mpiGrid, 1,+1,popID); + // update_remote_mapping_contribution(mpiGrid, 1,-1,popID); + //} else { + // update_remote_mapping_contribution_amr(mpiGrid, 1,+1,popID); + // update_remote_mapping_contribution_amr(mpiGrid, 1,-1,popID); + //} updateRemoteTimer.stop(); } @@ -316,7 +316,7 @@ void calculateSpatialTranslation( for (size_t c=0; cparameters[CellParams::LBWEIGHTCOUNTER] += time / localCells.size(); for (uint popID=0; popIDparameters[CellParams::LBWEIGHTCOUNTER] += mpiGrid[localCells[c]]->get_number_of_velocity_blocks(popID); + mpiGrid[localCells[c]]->parameters[CellParams::LBWEIGHTCOUNTER] = 1; // += mpiGrid[localCells[c]]->get_number_of_velocity_blocks(popID); } } } else {