Vamr final

datareduction/datareductionoperator.cpp

@@ -573,6 +573,9 @@ namespace DRO {
    // Blocks
    Blocks::Blocks(cuint _popID): DataReductionOperator(),popID(_popID) {
       popName=getObjectWrapper().particleSpecies[popID].name;
+      if (P::activateVamr) {
+        popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
    }
    Blocks::~Blocks() { }
 
@@ -633,6 +636,9 @@ namespace DRO {
    // and VariablePTensorOffDiagonal (23, 13, 12)
    VariablePTensorDiagonal::VariablePTensorDiagonal(cuint _popID): DataReductionOperator(),popID(_popID) {
       popName = getObjectWrapper().particleSpecies[popID].name;
+      if (P::activateVamr) {
+        popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
    }
    VariablePTensorDiagonal::~VariablePTensorDiagonal() { }
 
@@ -697,6 +703,9 @@ namespace DRO {
 
    VariablePTensorOffDiagonal::VariablePTensorOffDiagonal(cuint _popID): DataReductionOperator(),popID(_popID) {
       popName = getObjectWrapper().particleSpecies[popID].name;
+      if (P::activateVamr) {
+        popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
    }
    VariablePTensorOffDiagonal::~VariablePTensorOffDiagonal() { }
 
@@ -1030,6 +1039,9 @@ namespace DRO {
    // Rho nonthermal:
    VariableRhoNonthermal::VariableRhoNonthermal(cuint _popID): DataReductionOperator(),popID(_popID) {
       popName = getObjectWrapper().particleSpecies[popID].name;
+      if (P::activateVamr) {
+	    popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
       doSkip = (getObjectWrapper().particleSpecies[popID].thermalRadius == 0.0) ? true : false;
    }
    VariableRhoNonthermal::~VariableRhoNonthermal() { }
@@ -1059,6 +1071,9 @@ namespace DRO {
    // Rho thermal:
    VariableRhoThermal::VariableRhoThermal(cuint _popID): DataReductionOperator(),popID(_popID) {
       popName = getObjectWrapper().particleSpecies[popID].name;
+      if (P::activateVamr) {
+	    popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
       doSkip = (getObjectWrapper().particleSpecies[popID].thermalRadius == 0.0) ? true : false;
    }
    VariableRhoThermal::~VariableRhoThermal() { }
@@ -1088,6 +1103,9 @@ namespace DRO {
    // v nonthermal:
    VariableVNonthermal::VariableVNonthermal(cuint _popID): DataReductionOperator(),popID(_popID) {
       popName = getObjectWrapper().particleSpecies[popID].name;
+	  if (P::activateVamr) {
+        popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
       doSkip = (getObjectWrapper().particleSpecies[popID].thermalRadius == 0.0) ? true : false;
    }
    VariableVNonthermal::~VariableVNonthermal() { }
@@ -1122,6 +1140,9 @@ namespace DRO {
    //v thermal:
    VariableVThermal::VariableVThermal(cuint _popID): DataReductionOperator(),popID(_popID) {
       popName = getObjectWrapper().particleSpecies[popID].name;
+      if (P::activateVamr) {
+        popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
       doSkip = (getObjectWrapper().particleSpecies[popID].thermalRadius == 0.0) ? true : false;
    }
    VariableVThermal::~VariableVThermal() { }
@@ -1161,6 +1182,9 @@ namespace DRO {
    // and VariablePTensorNonthermalOffDiagonal (23, 13, 12)
    VariablePTensorNonthermalDiagonal::VariablePTensorNonthermalDiagonal(cuint _popID): DataReductionOperator(),popID(_popID) {
       popName = getObjectWrapper().particleSpecies[popID].name;
+      if (P::activateVamr) {
+        popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
       doSkip = (getObjectWrapper().particleSpecies[popID].thermalRadius == 0.0) ? true : false;
    }
    VariablePTensorNonthermalDiagonal::~VariablePTensorNonthermalDiagonal() { }
@@ -1209,6 +1233,9 @@ namespace DRO {
    // and VariablePTensorThermalOffDiagonal (23, 13, 12)
    VariablePTensorThermalDiagonal::VariablePTensorThermalDiagonal(cuint _popID): DataReductionOperator(),popID(_popID) {
       popName = getObjectWrapper().particleSpecies[popID].name;
+      if (P::activateVamr) {
+        popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
       doSkip = (getObjectWrapper().particleSpecies[popID].thermalRadius == 0.0) ? true : false;
    }
    VariablePTensorThermalDiagonal::~VariablePTensorThermalDiagonal() { }
@@ -1251,6 +1278,9 @@ namespace DRO {
 
    VariablePTensorNonthermalOffDiagonal::VariablePTensorNonthermalOffDiagonal(cuint _popID): DataReductionOperator(),popID(_popID) {
       popName = getObjectWrapper().particleSpecies[popID].name;
+	  if (P::activateVamr) {
+        popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
       doSkip = (getObjectWrapper().particleSpecies[popID].thermalRadius == 0.0) ? true : false;
    }
    VariablePTensorNonthermalOffDiagonal::~VariablePTensorNonthermalOffDiagonal() { }
@@ -1295,6 +1325,9 @@ namespace DRO {
 
    VariablePTensorThermalOffDiagonal::VariablePTensorThermalOffDiagonal(cuint _popID): DataReductionOperator(),popID(_popID) {
       popName = getObjectWrapper().particleSpecies[popID].name;
+      if (P::activateVamr) {
+        popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
       doSkip = (getObjectWrapper().particleSpecies[popID].thermalRadius == 0.0) ? true : false;
    }
    VariablePTensorThermalOffDiagonal::~VariablePTensorThermalOffDiagonal() { }
@@ -1340,6 +1373,9 @@ namespace DRO {
 
    VariableEffectiveSparsityThreshold::VariableEffectiveSparsityThreshold(cuint _popID): DataReductionOperator(),popID(_popID) {
      popName=getObjectWrapper().particleSpecies[popID].name;
+     if (P::activateVamr) {
+       popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+     }
    }
    VariableEffectiveSparsityThreshold::~VariableEffectiveSparsityThreshold() { }
 
@@ -1380,6 +1416,9 @@ namespace DRO {
     */
    VariablePrecipitationDiffFlux::VariablePrecipitationDiffFlux(cuint _popID): DataReductionOperatorHasParameters(),popID(_popID) {
       popName = getObjectWrapper().particleSpecies[popID].name;
+      if (P::activateVamr) {
+        popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
       lossConeAngle = getObjectWrapper().particleSpecies[popID].precipitationLossConeAngle; // deg
       emin = getObjectWrapper().particleSpecies[popID].precipitationEmin;    // already converted to SI
       emax = getObjectWrapper().particleSpecies[popID].precipitationEmax;    // already converted to SI
@@ -1505,6 +1544,9 @@ namespace DRO {
     */
    VariableMuSpace::VariableMuSpace(cuint _popID): DataReductionOperatorHasParameters(),popID(_popID) {
       popName = getObjectWrapper().particleSpecies[popID].name;
+      if (P::activateVamr) {
+        popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
       nBins = Parameters::PADmubins; //Number of bins to build muSpace
    }
    VariableMuSpace::~VariableMuSpace() { }
@@ -1594,6 +1636,9 @@ namespace DRO {
     */
    VariablePrecipitationLineDiffFlux::VariablePrecipitationLineDiffFlux(cuint _popID): DataReductionOperatorHasParameters(),popID(_popID) {
       popName = getObjectWrapper().particleSpecies[popID].name;
+      if (P::activateVamr) {
+	    popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
       emin = getObjectWrapper().particleSpecies[popID].precipitationEmin;    // already converted to SI
       emax = getObjectWrapper().particleSpecies[popID].precipitationEmax;    // already converted to SI
       nChannels = getObjectWrapper().particleSpecies[popID].precipitationNChannels; // number of energy channels, logarithmically spaced between emin and emax
@@ -1737,6 +1782,9 @@ namespace DRO {
     */
    VariableEnergyDensity::VariableEnergyDensity(cuint _popID): DataReductionOperatorHasParameters(),popID(_popID) {
       popName = getObjectWrapper().particleSpecies[popID].name;
+      if (P::activateVamr) {
+	    popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
       // Store internally in SI units
       solarwindenergy = getObjectWrapper().particleSpecies[popID].SolarWindEnergy;
       E1limit = solarwindenergy * getObjectWrapper().particleSpecies[popID].EnergyDensityLimit1;
@@ -1828,6 +1876,9 @@ namespace DRO {
    // q_i = m/2 * integral((v - <V>)^2 (v - <V>)_i * f(r,v) dV)
    VariableHeatFluxVector::VariableHeatFluxVector(cuint _popID): DataReductionOperator(),popID(_popID) {
       popName = getObjectWrapper().particleSpecies[popID].name;
+      if (P::activateVamr) {
+        popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
    }
    VariableHeatFluxVector::~VariableHeatFluxVector() { }
 
@@ -1913,6 +1964,9 @@ namespace DRO {
     */
    VariableNonMaxwellianity::VariableNonMaxwellianity(cuint _popID) : DataReductionOperator(), popID(_popID) {
       popName = getObjectWrapper().particleSpecies[popID].name;
+      if (P::activateVamr) {
+        popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
    }
    VariableNonMaxwellianity::~VariableNonMaxwellianity() {}
 

grid.cpp

@@ -285,8 +285,13 @@ void initializeGrids(
          mpiGrid[cells[i]]->parameters[CellParams::LBWEIGHTCOUNTER] = 0;
       }
 
+      if (P::activateVamr){
+		 //Fixed ghost to 1 for all the initial velocity cells to avoid destruction
+	     FixedGhost(mpiGrid,cells);
+      } 
+
       for (uint popID = 0; popID < getObjectWrapper().particleSpecies.size(); ++popID) {
-         adjustVelocityBlocks(mpiGrid, cells, true, popID);
+	     adjustVelocityBlocks(mpiGrid,cells,true,popID,true);
          // set initial LB metric based on number of blocks
          #pragma omp parallel for schedule(static)
          for (size_t i = 0; i < cells.size(); ++i) {
@@ -830,7 +835,8 @@ void prepareAMRLists(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGr
 bool adjustVelocityBlocks(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid,
                           const vector<CellID>& cellsToAdjust,
                           bool doPrepareToReceiveBlocks,
-                          const uint popID) {
+                          const uint popID,
+			              bool useGhost) {
    phiprof::Timer readjustBlocksTimer {"re-adjust blocks", {"Block adjustment"}};
    SpatialCell::setCommunicatedSpecies(popID);
 
@@ -844,7 +850,7 @@ bool adjustVelocityBlocks(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>&
    }
 
    // Batch call
-   update_velocity_block_content_lists(mpiGrid,validCells, popID);
+   update_velocity_block_content_lists(mpiGrid,validCells,useGhost,popID);
 
    // Get updated lists for blocks with content in spatial neighbours
    phiprof::Timer transferTimer {"Transfer with_content_list", {"MPI"}};

grid.h

@@ -135,7 +135,8 @@ void deallocateRemoteCellBlocks(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geomet
 bool adjustVelocityBlocks(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
                           const std::vector<CellID>& cellsToAdjust,
                           bool doPrepareToReceiveBlocks,
-                          const uint popID);
+                          const uint popID,
+			              bool useGhost);
 
 /*! Shrink to fit velocity space data to save memory.
  * \param mpiGrid Spatial grid

ioread.cpp

@@ -387,7 +387,10 @@ bool _readBlockDataCompressionNone(vlsv::ParallelReader & file,
    uint64_t byteSize;
    list<pair<string, string>> avgAttribs;
    bool success = true;
-   const string popName = getObjectWrapper().particleSpecies[popID].name;
+   string popName = getObjectWrapper().particleSpecies[popID].name;
+   if (P::activateVamr) {
+     popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+   }
    const string tagName = "BLOCKIDS";
 
    avgAttribs.push_back(make_pair("mesh", spatMeshName));
@@ -1028,7 +1031,10 @@ bool readBlockData(vlsv::ParallelReader& file, const string& meshName, const vec
    uint64_t* offsetArray = new uint64_t[N_processes];
 
    for (uint popID = 0; popID < getObjectWrapper().particleSpecies.size(); ++popID) {
-      const string& popName = getObjectWrapper().particleSpecies[popID].name;
+      string popName = getObjectWrapper().particleSpecies[popID].name;
+      if (P::activateVamr) {
+	     popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+      }
 
       // Create a cellID remapping lambda that can renumber our velocity space, should its size have changed.
       // By default, this is a no-op that keeps the blockIDs untouched.

iowrite.cpp

@@ -165,7 +165,10 @@ bool writeVelocityDistributionData(const uint popID, Writer& vlsvWriter,
    // In restart we just write velocity grids for all cells.
    // First write global Ids of those cells which write velocity blocks (here: all cells):
    map<string, string> attribs;
-   const string popName = getObjectWrapper().particleSpecies[popID].name;
+   string popName      = getObjectWrapper().particleSpecies[popID].name;
+   if (P::activateVamr) {
+     popName += std::to_string(getObjectWrapper().particleSpecies[popID].RefinementLevel);
+   }
    const string spatMeshName = "SpatialGrid";
    attribs["name"] = popName;
    bool success = true;
@@ -177,7 +180,6 @@ bool writeVelocityDistributionData(const uint popID, Writer& vlsvWriter,
       totalBlocks += mpiGrid[cells[i]]->get_number_of_velocity_blocks(popID);
       blocksPerCell.push_back(mpiGrid[cells[i]]->get_number_of_velocity_blocks(popID));
    }
-
    // The name of the mesh is "SpatialGrid"
    attribs["mesh"] = spatMeshName;
 
@@ -203,7 +205,7 @@ bool writeVelocityDistributionData(const uint popID, Writer& vlsvWriter,
    bbox[5] = vmesh::getMeshWrapper()->velocityMeshes->at(meshID).blockLength[2];
 
    attribs.clear();
-   attribs["mesh"] = getObjectWrapper().particleSpecies[popID].name;
+   attribs["mesh"] = popName;
    attribs["type"] = vlsv::mesh::STRING_UCD_AMR;
 
    // stringstream is necessary here to correctly convert refLevelMaxAllowed (hardcoded to zero now) into a string

object_wrapper.cpp

@@ -65,7 +65,10 @@ bool ObjectWrapper::addPopulationParameters() {
      RP::add(pop + "_vspace.vx_length","Initial number of velocity blocks in vx-direction.",1);
      RP::add(pop + "_vspace.vy_length","Initial number of velocity blocks in vy-direction.",1);
      RP::add(pop + "_vspace.vz_length","Initial number of velocity blocks in vz-direction.",1);
-     RP::add(pop + "_vspace.max_refinement_level","Maximum allowed mesh refinement level.", 1);
+     RP::add(pop + "_vspace.max_refinement_level","Old maximum allowed mesh refinement level.", 1);
+     RP::add(pop + "_vspace.vamr_refinement_level","New maximum allowed mesh refinement level. 0 is an homogenous grid, X will give X+1 grids.", 0);
+     RP::add(pop + "_vspace.vamr_criteria_method","Choice of the method for the vamr criteria. 0 : d><population>_sparse.minValue ,  1 : d > eps  , 2 : d > eps 2^-R . Eitheir 0 will be taken.", 0);
+     RP::add(pop + "_vspace.vamr_criteria_value","Value of epsilon (eps) for the vamr refinement criteria. Used only if vamr_criteria_method is 1 or 2.", 1e-15);
 
      // Thermal / suprathermal parameters
      Readparameters::add(pop + "_thermal.vx", "Center coordinate for the maxwellian distribution. Used for calculating the suprathermal moments.", -500000.0);
@@ -93,6 +96,44 @@ bool ObjectWrapper::addPopulationParameters() {
 bool ObjectWrapper::getPopulationParameters() {
    typedef Readparameters RP;
 
+   //Creation of the new species for each vAMR level
+   if(P::activateVamr) {
+     const int nbpopinit = getObjectWrapper().particleSpecies.size();
+     int shift=0;
+
+     for(int popID=0; popID < nbpopinit; popID++) {
+       species::Species& species=getObjectWrapper().particleSpecies[popID+shift];
+       vmesh::MeshParameters& vMesh=vmesh::getMeshWrapper()->velocityMeshesCreation->at(popID+shift);      
+       const std::string& pop = species.name;
+       RP::get(pop + "_vspace.vamr_refinement_level",  species.MaxRefinementLevel);
+       species.velocityMesh=popID+shift;
+       if ( species.MaxRefinementLevel > 0) {
+	     species.RefinementLevel=0;
+	     for (int R=0; R< species.MaxRefinementLevel; ++R) {
+	       species::Species newSpecies;
+	       vmesh::MeshParameters newVMesh;
+
+	   	   newSpecies.name = newVMesh.name = pop; // + std::to_string(R+1); Will only be done for the output files
+	   	   newSpecies.velocityMesh = popID+shift+1;
+		   //Insertion of the new level in order to have continuity for the same population in the future popID loop
+	   	   getObjectWrapper().particleSpecies.insert(getObjectWrapper().particleSpecies.begin()+popID+shift+1, newSpecies);
+	   	   vmesh::getMeshWrapper()->velocityMeshesCreation->insert(vmesh::getMeshWrapper()->velocityMeshesCreation->begin() +popID+shift+1,newVMesh);
+
+	   	   species::Species& species2=getObjectWrapper().particleSpecies[popID+shift+1];
+	   	   vmesh::MeshParameters& vMesh2=vmesh::getMeshWrapper()->velocityMeshesCreation->at(popID+shift+1);
+
+	   	   species2.RefinementLevel=R+1;
+	   	   species2.MaxRefinementLevel= species.MaxRefinementLevel;
+
+	   	   shift+=1;
+	     }
+       }else {
+	     species.MaxRefinementLevel=0;
+	     species.RefinementLevel=0;
+       }
+     }
+   }
+
    // Particle population parameters
    for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) {
 
@@ -144,6 +185,20 @@ bool ObjectWrapper::getPopulationParameters() {
       RP::get(pop + "_vspace.vx_length",vMesh.gridLength[0]);
       RP::get(pop + "_vspace.vy_length",vMesh.gridLength[1]);
       RP::get(pop + "_vspace.vz_length",vMesh.gridLength[2]);
+
+      // Vamr parameters that are saved in species
+      RP::get(pop + "_vspace.vamr_criteria_method",  species.CriteriaMethod);
+      RP::get(pop + "_vspace.vamr_criteria_value",  species.CriteriaValue);
+
+      if(P::activateVamr && species.RefinementLevel>0) {
+		//Nv(R)=Nv(R=0)*2^R
+	    vMesh.gridLength[0] *= (1u << species.RefinementLevel);
+	    vMesh.gridLength[1] *= (1u << species.RefinementLevel);
+	    vMesh.gridLength[2] *= (1u << species.RefinementLevel);
+	    //species.sparseMinValue *= std::pow(species.minValueRefinementShift, species.RefinementLevel);
+        species.sparseBlockAddWidthV=0;
+      }
+
       if(vMesh.gridLength[0] > MAX_BLOCKS_PER_DIM  ||
             vMesh.gridLength[1] > MAX_BLOCKS_PER_DIM  ||
             vMesh.gridLength[2] > MAX_BLOCKS_PER_DIM ) {

parameters.cpp

@@ -129,6 +129,9 @@ bool P::propagateVlasovTranslation = true;
 bool P::propagateField = true;
 
 bool P::dynamicTimestep = true;
+bool P::activateVamr = false;
+uint P::vAMRorder = 1;
+uint P::vAMRrefineStep = 5;
 
 Real P::maxWaveVelocity = 0.0;
 uint P::maxFieldSolverSubcycles = 0.0;
@@ -325,6 +328,9 @@ bool P::addParameters() {
            "zero length timesteps.",
            true);
    RP::add("dynamic_timestep", "If true,  timestep is set based on  CFL limits (default on)", true);
+   RP::add("activate_vamr","Activate the velocity mesh refinement. If 1 it is activated. ", false);
+   RP::add("vamr_refineStep","Number of steps between each velocity refinement check.", 5);
+   RP::add("vamr_order","Select the numerical order for the creation of cells in the vamr. Existing orders :1, 3 and 5 , either 1 will be taken.", 1);
    RP::add("hallMinimumRho",
            "Minimum rho value used for the Hall and electron pressure gradient terms in the Lorentz force and in the "
            "field solver. Default is very low and has no effect in practice.",
@@ -825,6 +831,9 @@ void Parameters::getParameters() {
    RP::get("propagate_vlasov_acceleration", P::propagateVlasovAcceleration);
    RP::get("propagate_vlasov_translation", P::propagateVlasovTranslation);
    RP::get("dynamic_timestep", P::dynamicTimestep);
+   RP::get("activate_vamr", P::activateVamr);
+   RP::get("vamr_order", P::vAMRorder);
+   RP::get("vamr_refineStep", P::vAMRrefineStep);
    Real hallRho;
    RP::get("hallMinimumRho", hallRho);
    P::hallMinimumRhom = hallRho * physicalconstants::MASS_PROTON;