Multiple AMR boxes

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

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"<<std::endl;
-      }
-   }
    globalflags::AMRstencilWidth = neighborhood_size;
 
    const std::array<uint64_t, 3> 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<SpatialCell,dccrg::Cartesian_Geometry>& 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"<<std::endl;
-      }
-   }
    int full_neighborhood_size = max(2, VLASOV_STENCIL_WIDTH);
 
    neighborhood.clear();
@@ -1039,15 +1003,6 @@ void initializeStencils(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& 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";

iowrite.cpp

@@ -1104,7 +1104,7 @@ bool writeVelocitySpace(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& 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;

memoryallocation.h

@@ -24,6 +24,7 @@
 
 #include <cstdlib>
 #include <cstddef>
+#include <cstdint>
 #include <stdexcept>
 #include <string.h>
 

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<uint> P::amrBoxHalfWidthX;
+std::vector<uint> P::amrBoxHalfWidthY;
+std::vector<uint> P::amrBoxHalfWidthZ;
+std::vector<Realf> P::amrBoxCenterX;
+std::vector<Realf> P::amrBoxCenterY;
+std::vector<Realf> P::amrBoxCenterZ;
+std::vector<int> P::amrBoxMaxLevel;
 vector<string> P::blurPassString;
 std::vector<int> 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;
          }

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<uint> amrBoxHalfWidthX;
+   static std::vector<uint> amrBoxHalfWidthY;
+   static std::vector<uint> amrBoxHalfWidthZ;
+   static std::vector<Realf> amrBoxCenterX;
+   static std::vector<Realf> amrBoxCenterY;
+   static std::vector<Realf> amrBoxCenterZ;
+   static std::vector<int> amrBoxMaxLevel;
    static bool amrTransShortPencils;        /*!< Use short or longpencils in AMR translation.*/
    static std::vector<std::string> blurPassString;
    static std::vector<int> numPasses;

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

projects/Flowthrough/Flowthrough.cpp

@@ -266,14 +266,16 @@ namespace projects {
       std::vector<CellID> cells = mpiGrid.get_cells();
       for (CellID id : cells) {
          std::array<double,3> 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);
+            }
          }
       }
 

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