diff --git a/datareduction/datareducer.cpp b/datareduction/datareducer.cpp index 0e3a300418..6ffd4b4265 100644 --- a/datareduction/datareducer.cpp +++ b/datareduction/datareducer.cpp @@ -239,7 +239,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti if(P::systemWriteAllDROs || lowercase == "populations_rho" || lowercase == "populations_vg_rho") { // Per-population particle number density for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::DataReductionOperatorPopulations(pop + "/vg_rho", i, offsetof(spatial_cell::Population, RHO), 1)); outputReducer->addMetadata(outputReducer->size()-1,"1/m^3","$\\mathrm{m}^{-3}$","$n_\\mathrm{"+pop+"}$","1.0"); } @@ -289,7 +292,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti if(P::systemWriteAllDROs || lowercase == "populations_v" || lowercase == "populations_vg_v") { // Per population bulk velocities for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::DataReductionOperatorPopulations(pop + "/vg_v", i, offsetof(spatial_cell::Population, V), 3)); outputReducer->addMetadata(outputReducer->size()-1,"m/s","$\\mathrm{m}\\,\\mathrm{s}^{-1}$","$V_\\mathrm{"+pop+"}$","1.0"); } @@ -300,7 +306,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti if(P::systemWriteAllDROs || lowercase == "populations_moments_backstream" || lowercase == "populations_moments_nonthermal" || lowercase == "populations_vg_moments_nonthermal") { // Per-population moments of the backstreaming part for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::VariableRhoNonthermal(i)); outputReducer->addOperator(new DRO::VariableVNonthermal(i)); outputReducer->addOperator(new DRO::VariablePTensorNonthermalDiagonal(i)); @@ -317,7 +326,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti if(P::systemWriteAllDROs || lowercase == "populations_moments_nonbackstream" || lowercase == "populations_moments_thermal" || lowercase == "populations_vg_moments_thermal") { // Per-population moments of the non-backstreaming (thermal?) part. for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::VariableRhoThermal(i)); outputReducer->addOperator(new DRO::VariableVThermal(i)); outputReducer->addOperator(new DRO::VariablePTensorThermalDiagonal(i)); @@ -335,7 +347,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti // Effective sparsity threshold affecting each cell, if dynamic threshould algorithm is used for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::VariableEffectiveSparsityThreshold(i)); outputReducer->addMetadata(outputReducer->size()-1,"s^3/m^6","$\\mathrm{m}^{-6}\\,\\mathrm{s}^{3}$","$f_\\mathrm{"+pop+",min}$","1.0"); } @@ -347,7 +362,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti // Accumulated lost particle number, per population, in each cell, since last restart for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::DataReductionOperatorPopulations(pop + "/vg_rho_loss_adjust", i, offsetof(spatial_cell::Population, RHOLOSSADJUST), 1)); outputReducer->addMetadata(outputReducer->size()-1,"1/m^3","$\\mathrm{m}^{-3}$","$\\Delta_\\mathrm{loss} n_\\mathrm{"+pop+"}$","1.0"); } @@ -375,7 +393,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti // Per-population maximum timestep constraint as calculated by the velocity space vlasov update for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::DataReductionOperatorPopulations(pop + "/vg_maxdt_acceleration", i, offsetof(spatial_cell::Population, max_dt[1]), 1)); outputReducer->addMetadata(outputReducer->size()-1,"s","$\\mathrm{s}$","$\\Delta t_\\mathrm{"+pop+",V,max}$","1.0"); } @@ -395,7 +416,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti // Per-population maximum timestep constraint as calculated by the real space vlasov update for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::DataReductionOperatorPopulations(pop + "/vg_maxdt_translation", i, offsetof(spatial_cell::Population, max_dt[0]), 1)); outputReducer->addMetadata(outputReducer->size()-1,"s","$\\mathrm{s}$","$\\Delta t_\\mathrm{"+pop+",R,max}$","1.0"); } @@ -407,7 +431,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti // Per-population energy density in three energy ranges for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::VariableEnergyDensity(i)); std::stringstream conversion; conversion << (1.0e-6)/physicalconstants::CHARGE; @@ -421,7 +448,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti // Per-population precipitation differential flux (within loss cone) for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::VariablePrecipitationDiffFlux(i)); std::stringstream conversion; conversion << (1.0e-4)*physicalconstants::CHARGE; @@ -435,7 +465,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti // Per-population 1d muspace for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::VariableMuSpace(i)); outputReducer->addMetadata(outputReducer->size()-1,"1/m^3","$\\mathrm{m}^{-3}$","$f(\\mu)_\\mathrm{"+pop+"}$","1.0"); } @@ -447,7 +480,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti // Per-population precipitation differential flux (along line) for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::VariablePrecipitationLineDiffFlux(i)); std::stringstream conversion; conversion << (1.0e-4)*physicalconstants::CHARGE; @@ -461,7 +497,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti // Per-population heat flux vector for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::VariableHeatFluxVector(i)); outputReducer->addMetadata(outputReducer->size()-1,"W/m^2","$\\mathrm{W}\\,\\mathrm{m}^{-2}$","$q_\\mathrm{"+pop+"}$","1.0"); } @@ -473,7 +512,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti // Per-population dimensionless non-maxwellianity parameter for (unsigned int i = 0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species = getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::VariableNonMaxwellianity(i)); outputReducer->addMetadata(outputReducer->size() - 1, "", "", "$\\tilde{\\epsilon}_\\mathrm{M," + pop + "}$", "1.0"); @@ -616,7 +658,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti // Per-population velocity space block counts for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::Blocks(i)); outputReducer->addMetadata(outputReducer->size()-1,"","","$\\mathrm{"+pop+" blocks}$",""); } @@ -660,7 +705,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti // Per-population number of subcycles performed for velocity space update for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::DataReductionOperatorPopulations(pop + "/vg_acceleration_subcycles", i, offsetof(spatial_cell::Population, ACCSUBCYCLES), 1)); outputReducer->addMetadata(outputReducer->size()-1,"","","$\\mathrm{"+pop+" Acc subcycles}$",""); } @@ -815,7 +863,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti // Per-population pressure tensor, stored as diagonal and offdiagonal components for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } outputReducer->addOperator(new DRO::VariablePTensorDiagonal(i)); outputReducer->addMetadata(outputReducer->size()-1,"Pa","$\\mathrm{Pa}$","$\\mathcal{P}_\\mathrm{"+pop+"}$","1.0"); outputReducer->addOperator(new DRO::VariablePTensorOffDiagonal(i)); @@ -2813,7 +2864,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti // Per-particle overall lost particle number for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } diagnosticReducer->addOperator(new DRO::DataReductionOperatorPopulations(pop + "/vg_rho_loss_adjust", i, offsetof(spatial_cell::Population, RHOLOSSADJUST), 1)); } if(!P::diagnosticWriteAllDROs) { @@ -2847,7 +2901,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti if(P::diagnosticWriteAllDROs || lowercase == "populations_maxrdt" || lowercase == "populations_maxdt_translation" || lowercase == "populations_vg_maxdt_translation") { for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } diagnosticReducer->addOperator(new DRO::DataReductionOperatorPopulations(pop + "/vg_maxdt_translation", i, offsetof(spatial_cell::Population, max_dt[0]), 1)); } if(!P::diagnosticWriteAllDROs) { @@ -2857,7 +2914,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti if(P::diagnosticWriteAllDROs || lowercase == "populations_maxvdt" || lowercase == "populations_maxdt_acceleration" || lowercase == "populations_vg_maxdt_acceleration") { for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) { species::Species& species=getObjectWrapper().particleSpecies[i]; - const std::string& pop = species.name; + std::string pop = species.name; + if (P::activateVamr) { + pop += std::to_string(species.RefinementLevel); + } diagnosticReducer->addOperator(new DRO::DataReductionOperatorPopulations(pop + "/vg_maxdt_acceleration", i, offsetof(spatial_cell::Population, max_dt[1]), 1)); } if(!P::diagnosticWriteAllDROs) { diff --git a/datareduction/datareductionoperator.cpp b/datareduction/datareductionoperator.cpp index 6b8f5c1150..4c3f37d0a0 100644 --- a/datareduction/datareductionoperator.cpp +++ b/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 - )^2 (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() {} diff --git a/grid.cpp b/grid.cpp index b82879cdf9..a89cba019a 100644 --- a/grid.cpp +++ b/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& mpiGr bool adjustVelocityBlocks(dccrg::Dccrg& mpiGrid, const vector& 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& } // 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"}}; diff --git a/grid.h b/grid.h index 32c5d297f5..0ce5ca3338 100644 --- a/grid.h +++ b/grid.h @@ -135,7 +135,8 @@ void deallocateRemoteCellBlocks(dccrg::Dccrg& mpiGrid, const std::vector& cellsToAdjust, bool doPrepareToReceiveBlocks, - const uint popID); + const uint popID, + bool useGhost); /*! Shrink to fit velocity space data to save memory. * \param mpiGrid Spatial grid diff --git a/ioread.cpp b/ioread.cpp index 0a72ce7661..c338325a79 100644 --- a/ioread.cpp +++ b/ioread.cpp @@ -387,7 +387,10 @@ bool _readBlockDataCompressionNone(vlsv::ParallelReader & file, uint64_t byteSize; list> 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. diff --git a/iowrite.cpp b/iowrite.cpp index d8aaa1112d..ee2d829d08 100644 --- a/iowrite.cpp +++ b/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 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 diff --git a/object_wrapper.cpp b/object_wrapper.cpp index a7a5afdf5e..abc4f4bf65 100644 --- a/object_wrapper.cpp +++ b/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>_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 ) { diff --git a/parameters.cpp b/parameters.cpp index 792710e76d..279e4e3ff8 100644 --- a/parameters.cpp +++ b/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; diff --git a/parameters.h b/parameters.h index 8f0c59cca5..4f24e73e9c 100644 --- a/parameters.h +++ b/parameters.h @@ -136,6 +136,9 @@ struct Parameters { static bool propagateVlasovTranslation; /*!< If true, distribution function is propagated in ordinary space during the simulation.*/ + static bool activateVamr; /*!< If true, the velocity mesh refinement is activated.*/ + static uint vAMRorder; /*!< Order of the vAMR method */ + static uint vAMRrefineStep; /*!< Number of step between two vAMR grid refinement*/ static Real maxWaveVelocity; /*!< Maximum wave velocity allowed in LDZ. */ static uint maxFieldSolverSubcycles; /*!< Maximum allowed field solver subcycles. */ static Real resistivity; /*!< Resistivity in Ohm's law eta*J term. */ diff --git a/particle_species.cpp b/particle_species.cpp index 4d71c32aa0..114e389d53 100644 --- a/particle_species.cpp +++ b/particle_species.cpp @@ -29,6 +29,6 @@ using namespace std; species::Species::Species() { } species::Species::Species(const Species& other) : - name {other.name}, charge {other.charge}, mass {other.mass}, sparseMinValue {other.sparseMinValue}, velocityMesh {other.velocityMesh} {} + name {other.name}, charge {other.charge}, mass {other.mass}, sparseMinValue {other.sparseMinValue}, velocityMesh {other.velocityMesh}, RefinementLevel {other.RefinementLevel}, MaxRefinementLevel {other.MaxRefinementLevel}, CriteriaMethod {other.CriteriaMethod}, CriteriaValue {other.CriteriaValue} {} species::Species::~Species() { } diff --git a/particle_species.h b/particle_species.h index 0bf4f2cd67..944056ebeb 100644 --- a/particle_species.h +++ b/particle_species.h @@ -68,6 +68,10 @@ namespace species { Real precipitationEmin; /*!< Lowest energy channel (in keV) for precipitation differential flux evaluation. Default 0.1. */ Real precipitationEmax; /*!< Highest energy channel (in keV) for precipitation differential flux evaluation. Default 100. */ Real precipitationLossConeAngle; /*!< Fixed loss cone opening angle (in deg) for precipitation differential flux evaluation. Default 10. */ + int RefinementLevel; /*!< The level of refinement of this vmesh with 0 the coarsest velocity grid */ + int MaxRefinementLevel; /*!< The maximum level of refinement of the velocity grid of this population */ + int CriteriaMethod; /*!< The method used for the vamr criteria */ + Real CriteriaValue; /*!< The value epsilon in the vamr criteria */ Species(); Species(const Species& other); diff --git a/spatial_cells/block_adjust_cpu.cpp b/spatial_cells/block_adjust_cpu.cpp index 86c91d8a78..3ccccafe19 100644 --- a/spatial_cells/block_adjust_cpu.cpp +++ b/spatial_cells/block_adjust_cpu.cpp @@ -33,6 +33,7 @@ namespace spatial_cell { void update_velocity_block_content_lists( dccrg::Dccrg& mpiGrid, const vector& cells, + bool useGhost, const uint popID) { if (cells.size()==0) { @@ -46,7 +47,7 @@ namespace spatial_cell { #pragma omp for schedule(dynamic) for (uint i=0; iupdateSparseMinValue(popID); - mpiGrid[cells[i]]->update_velocity_block_content_lists(popID); + mpiGrid[cells[i]]->update_velocity_block_content_lists(popID,useGhost); } // timer.stop(); } // end parallel region diff --git a/spatial_cells/block_adjust_cpu.hpp b/spatial_cells/block_adjust_cpu.hpp index 87357b03c6..88bd7e95bb 100644 --- a/spatial_cells/block_adjust_cpu.hpp +++ b/spatial_cells/block_adjust_cpu.hpp @@ -46,6 +46,7 @@ namespace spatial_cell { void update_velocity_block_content_lists( dccrg::Dccrg& mpiGrid, const vector& cells, + bool useGhost, const uint popID=0); void adjust_velocity_blocks_in_cells( diff --git a/spatial_cells/spatial_cell_cpu.cpp b/spatial_cells/spatial_cell_cpu.cpp index aab9ceaa7f..324f237bf4 100644 --- a/spatial_cells/spatial_cell_cpu.cpp +++ b/spatial_cells/spatial_cell_cpu.cpp @@ -277,7 +277,7 @@ namespace spatial_cell { //neighbor_ptrs is empty, so we adjust only based on local velocity space. std::vector neighbor_ptrs; - update_velocity_block_content_lists(popID); + update_velocity_block_content_lists(popID,false); adjust_velocity_blocks(neighbor_ptrs,popID,doDeleteEmpty); } @@ -287,7 +287,7 @@ namespace spatial_cell { sense in given block. Also returns false if given block doesn't exist or is an error block. */ - bool SpatialCell::compute_block_has_content(const vmesh::LocalID& blockLID,const uint popID) const { + bool SpatialCell::compute_block_has_content(const vmesh::LocalID& blockLID,const uint popID,bool useGhost) const { debug_population_check(popID); #ifdef DEBUG_SPATIAL_CELL const vmesh::GlobalID blockGID = populations[popID].vmesh->getGlobalID(blockLID); @@ -303,14 +303,36 @@ namespace spatial_cell { } bool has_content = false; - const Real velocity_block_min_value = getVelocityBlockMinValue(popID); - const Realf* block_data = populations[popID].blockContainer->getData(blockLID); - for (unsigned int i=0; i= velocity_block_min_value) { - has_content = true; - break; - } - } + + if (!P::activateVamr || getObjectWrapper().particleSpecies[popID].RefinementLevel==0 || !useGhost ){ + const Real velocity_block_min_value = getVelocityBlockMinValue(popID); + const Realf* block_data = populations[popID].blockContainer->getData(blockLID); + for (unsigned int i=0; i= velocity_block_min_value) { + has_content = true; + break; + } + } + }else{ + //Newly created cells that don't respect the vAMR criterion should be destroyed + uint8_t *ghost = populations[popID].blockContainer->getGhost(blockLID); + bool has_content1 = false; + if (ghost[0]==1){// ghost=2 may be used later for border + has_content1 = true; + } + bool has_content2 = false; + const Real velocity_block_min_value = getVelocityBlockMinValue(popID); + const Realf* block_data = populations[popID].blockContainer->getData(blockLID); + for (unsigned int i=0; i= velocity_block_min_value) { + has_content2 = true; + break; + } + } + if (has_content1 && has_content2){ + has_content = true; + } + } return has_content; } @@ -412,6 +434,12 @@ namespace spatial_cell { block_lengths.push_back(sizeof(Realf) * WID3 * populations[activePopID].blockContainer->size()); } + if ((SpatialCell::mpi_transfer_type & Transfer::VEL_BLOCK_REFINED) !=0) { + //Not used but can be to share the Refined parameter + displacements.push_back((uint8_t*) get_refined(activePopID) - (uint8_t*) this); + block_lengths.push_back(sizeof(uint8_t) * WID3 * populations[activePopID].blockContainer->size()); + } + if ((SpatialCell::mpi_transfer_type & Transfer::NEIGHBOR_VEL_BLOCK_DATA) != 0) { /*We are actually transferring the data of a * neighbor. The values of neighbor_block_data @@ -672,14 +700,14 @@ namespace spatial_cell { /** Update the two lists containing blocks with content, and blocks without content. * @see adjustVelocityBlocks */ - void SpatialCell::update_velocity_block_content_lists(const uint popID) { + void SpatialCell::update_velocity_block_content_lists(const uint popID,bool useGhost) { debug_population_check(popID); velocity_block_with_content_list->clear(); velocity_block_with_no_content_list->clear(); for (vmesh::LocalID block_index=0; block_indexsize(); ++block_index) { const vmesh::GlobalID globalID = populations[popID].vmesh->getGlobalID(block_index); - if (compute_block_has_content(block_index,popID)){ + if (compute_block_has_content(block_index,popID,useGhost)){ velocity_block_with_content_list->push_back(globalID); } else { velocity_block_with_no_content_list->push_back(globalID); diff --git a/spatial_cells/spatial_cell_cpu.hpp b/spatial_cells/spatial_cell_cpu.hpp index 65b72beee2..77ce548a90 100644 --- a/spatial_cells/spatial_cell_cpu.hpp +++ b/spatial_cells/spatial_cell_cpu.hpp @@ -251,6 +251,14 @@ namespace spatial_cell { const Realf* get_data(const uint popID) const; Realf* get_data(const vmesh::LocalID& blockLID,const uint popID); const Realf* get_data(const vmesh::LocalID& blockLID,const uint popID) const; + uint8_t* get_refined(const uint popID); + const uint8_t* get_refined(const uint popID) const; + uint8_t* get_refined(const vmesh::LocalID& blockLID,const uint popID); + const uint8_t* get_refined(const vmesh::LocalID& blockLID,const uint popID) const; + uint8_t* get_ghost(const uint popID); + const uint8_t* get_ghost(const uint popID) const; + uint8_t* get_ghost(const vmesh::LocalID& blockLID,const uint popID); + const uint8_t* get_ghost(const vmesh::LocalID& blockLID,const uint popID) const; Real* get_block_parameters(const uint popID); const Real* get_block_parameters(const uint popID) const; Real* get_block_parameters(const vmesh::LocalID& blockLID,const uint popID); @@ -324,7 +332,7 @@ namespace spatial_cell { // Templated function for storing a v-space read from a file or generated elsewhere template void add_velocity_blocks(const uint popID,const std::vector& blocks,fileReal* avgBuffer); - void update_velocity_block_content_lists(const uint popID); + void update_velocity_block_content_lists(const uint popID,bool useGhost); bool checkMesh(const uint popID); void clear(const uint popID, bool shrink=false); void setNewSizeClear(const uint popID, const vmesh::LocalID& newSize); @@ -385,7 +393,7 @@ namespace spatial_cell { private: //SpatialCell& operator=(const SpatialCell&); - bool compute_block_has_content(const vmesh::GlobalID& block,const uint popID) const; + bool compute_block_has_content(const vmesh::GlobalID& block,const uint popID,const bool useGhost) const; static int activePopID; bool initialized; @@ -445,6 +453,46 @@ namespace spatial_cell { return populations[popID].blockContainer->getData(blockLID); } + inline uint8_t* SpatialCell::get_refined(const uint popID) { + debug_population_check(popID); + return populations[popID].blockContainer->getRefined(); + } + + inline const uint8_t* SpatialCell::get_refined(const uint popID) const { + debug_population_check(popID); + return populations[popID].blockContainer->getRefined(); + } + + inline uint8_t* SpatialCell::get_refined(const vmesh::LocalID& blockLID,const uint popID) { + debug_population_check(popID,blockLID); + return populations[popID].blockContainer->getRefined(blockLID); + } + + inline const uint8_t* SpatialCell::get_refined(const vmesh::LocalID& blockLID,const uint popID) const { + debug_population_check(popID,blockLID); + return populations[popID].blockContainer->getRefined(blockLID); + } + + inline uint8_t* SpatialCell::get_ghost(const uint popID) { + debug_population_check(popID); + return populations[popID].blockContainer->getGhost(); + } + + inline const uint8_t* SpatialCell::get_ghost(const uint popID) const { + debug_population_check(popID); + return populations[popID].blockContainer->getGhost(); + } + + inline uint8_t* SpatialCell::get_ghost(const vmesh::LocalID& blockLID,const uint popID) { + debug_population_check(popID,blockLID); + return populations[popID].blockContainer->getGhost(blockLID); + } + + inline const uint8_t* SpatialCell::get_ghost(const vmesh::LocalID& blockLID,const uint popID) const { + debug_population_check(popID,blockLID); + return populations[popID].blockContainer->getGhost(blockLID); + } + inline Real* SpatialCell::get_block_parameters(const uint popID) { debug_population_check(popID); return populations[popID].blockContainer->getParameters(); diff --git a/spatial_cells/spatial_cell_wrapper.hpp b/spatial_cells/spatial_cell_wrapper.hpp index 65b9fb4228..7e647110dd 100644 --- a/spatial_cells/spatial_cell_wrapper.hpp +++ b/spatial_cells/spatial_cell_wrapper.hpp @@ -54,6 +54,7 @@ namespace spatial_cell { static const uint64_t VEL_BLOCK_LIST_STAGE1 = (1ull<<2); static const uint64_t VEL_BLOCK_LIST_STAGE2 = (1ull<<3); static const uint64_t VEL_BLOCK_DATA = (1ull<<4); + static const uint64_t VEL_BLOCK_REFINED = (1ull<<5);//Not used but can be to share the Refined parameter static const uint64_t VEL_BLOCK_PARAMETERS = (1ull<<6); static const uint64_t VEL_BLOCK_WITH_CONTENT_STAGE1 = (1ull<<7); static const uint64_t VEL_BLOCK_WITH_CONTENT_STAGE2 = (1ull<<8); diff --git a/spatial_cells/velocity_block_container.h b/spatial_cells/velocity_block_container.h index 0cdf58d3cb..e07a8fc824 100644 --- a/spatial_cells/velocity_block_container.h +++ b/spatial_cells/velocity_block_container.h @@ -73,6 +73,14 @@ namespace vmesh { ARCH_HOSTDEV const Realf* getData() const; ARCH_HOSTDEV Realf* getData(const vmesh::LocalID blockLID); ARCH_HOSTDEV const Realf* getData(const vmesh::LocalID blockLID) const; + ARCH_HOSTDEV uint8_t* getRefined(); + ARCH_HOSTDEV const uint8_t* getRefined() const; + ARCH_HOSTDEV uint8_t* getRefined(const vmesh::LocalID blockLID); + ARCH_HOSTDEV const uint8_t* getRefined(const vmesh::LocalID blockLID) const; + ARCH_HOSTDEV uint8_t* getGhost(); + ARCH_HOSTDEV const uint8_t* getGhost() const; + ARCH_HOSTDEV uint8_t* getGhost(const vmesh::LocalID blockLID); + ARCH_HOSTDEV const uint8_t* getGhost(const vmesh::LocalID blockLID) const; ARCH_HOSTDEV Real* getParameters(); ARCH_HOSTDEV const Real* getParameters() const; ARCH_HOSTDEV Real* getParameters(const vmesh::LocalID blockLID); @@ -104,8 +112,12 @@ namespace vmesh { #ifdef DEBUG_VBC const Realf& getData(const vmesh::LocalID blockLID,const unsigned int cell) const; + const uint8_t& getRefined(const vmesh::LocalID blockLID,const unsigned int cell) const; + const uint8_t& getGhost(const vmesh::LocalID blockLID,const unsigned int cell) const; const Real& getParameters(const vmesh::LocalID blockLID,const unsigned int i) const; void setData(const vmesh::LocalID blockLID,const unsigned int cell,const Realf value); + void setRefined(const vmesh::LocalID blockLID,const unsigned int cell,const uint8_t value); + void setGhost(const vmesh::LocalID blockLID,const unsigned int cell,const uint8_t value); #endif private: @@ -114,11 +126,15 @@ namespace vmesh { #ifdef USE_GPU split::SplitVector block_data; + split::SplitVector refined; + split::SplitVector ghost; split::SplitVector parameters; size_t cachedCapacity; size_t cachedSize; #else std::vector > block_data; + std::vector refined; + std::vector ghost; std::vector > parameters; #endif }; @@ -126,15 +142,21 @@ namespace vmesh { inline VelocityBlockContainer::VelocityBlockContainer() { #ifdef USE_GPU block_data = split::SplitVector(INIT_VMESH_SIZE*WID3); + refined = split::SplitVector(INIT_VMESH_SIZE*WID3); + ghost = split::SplitVector(INIT_VMESH_SIZE); parameters = split::SplitVector(INIT_VMESH_SIZE*BlockParams::N_VELOCITY_BLOCK_PARAMS); cachedCapacity = INIT_VMESH_SIZE; cachedSize = 0; #else block_data = std::vector>(WID3); + refined = std::vector(WID3); + ghost = std::vector(); parameters = std::vector>(BlockParams::N_VELOCITY_BLOCK_PARAMS); //cachedCapacity = 1; #endif block_data.clear(); + refined.clear(); + ghost.clear(); parameters.clear(); // gpuStream_t stream = gpu_getStream(); } @@ -144,15 +166,21 @@ namespace vmesh { inline VelocityBlockContainer::VelocityBlockContainer(const VelocityBlockContainer& other) { #ifdef USE_GPU block_data = split::SplitVector(other.cachedCapacity*WID3); + refined = split::SplitVector(other.cachedCapacity*WID3); + ghost = split::SplitVector(other.cachedCapacity); parameters = split::SplitVector(other.cachedCapacity*BlockParams::N_VELOCITY_BLOCK_PARAMS); // Overwrite is like a copy assign but takes a stream gpuStream_t stream = gpu_getStream(); block_data.overwrite(other.block_data,stream); + refined.overwrite(other.refined,stream); + ghost.overwrite(other.ghost,stream); parameters.overwrite(other.parameters,stream); cachedSize = other.cachedSize; cachedCapacity = other.cachedCapacity; #else block_data = std::vector>(other.block_data); + refined = std::vector(other.refined); + ghost = std::vector(other.ghost); parameters = std::vector>(other.parameters); // block_data.reserve(other.capacity()*WID3); // parameters.reserve(other.capacity()*BlockParams::N_VELOCITY_BLOCK_PARAMS); @@ -163,14 +191,20 @@ namespace vmesh { #ifdef USE_GPU gpuStream_t stream = gpu_getStream(); block_data.reserve(other.cachedCapacity*WID3, true, stream); + refined.reserve(other.cachedCapacity*WID3, true, stream); + ghost.reserve(other.cachedCapacity, true, stream); parameters.reserve(other.cachedCapacity*BlockParams::N_VELOCITY_BLOCK_PARAMS, true, stream); // Overwrite is like a copy assign but takes a stream block_data.overwrite(other.block_data,stream); + refined.overwrite(other.refined,stream); + ghost.overwrite(other.ghost,stream); parameters.overwrite(other.parameters,stream); cachedSize = other.cachedSize; cachedCapacity = other.cachedCapacity; #else block_data = other.block_data; + refined = other.refined; + ghost = other.ghost; parameters = other.parameters; // block_data.reserve(other.capacity()*WID3); // parameters.reserve(other.capacity()*BlockParams::N_VELOCITY_BLOCK_PARAMS); @@ -222,18 +256,24 @@ namespace vmesh { if (shrink) { cachedCapacity = 1; block_data = split::SplitVector(WID3); + refined = split::SplitVector(WID3); + ghost = split::SplitVector(); parameters = split::SplitVector(BlockParams::N_VELOCITY_BLOCK_PARAMS); } #else if (shrink) { block_data = std::vector>(WID3); + refined = std::vector(WID3); + ghost = std::vector(); parameters = std::vector>(BlockParams::N_VELOCITY_BLOCK_PARAMS); } #endif block_data.clear(); + ghost.clear(); + refined.clear(); parameters.clear(); #ifdef DEBUG_VBC - if ((block_data.size() != 0) || (parameters.size() != 0)) { + if ((block_data.size() != 0) || (refined.size() != 0) || (ghost.size() != 0) || (parameters.size() != 0)) { std::cerr<<"VBC CLEAR FAILED"<= numberOfBlocks) ok = false; if (source >= currentCapacity) ok = false; + if (source >= numberOfBlocksR) ok = false; + if (source >= currentCapacityR) ok = false; + if (source >= numberOfBlocksG) ok = false; + if (source >= currentCapacityG) ok = false; if (source >= numberOfBlocksP) ok = false; if (source >= currentCapacityP) ok = false; if (target >= numberOfBlocks) ok = false; @@ -265,22 +313,27 @@ namespace vmesh { if (numberOfBlocks > currentCapacity) ok = false; if (source != numberOfBlocks-1) ok = false; // only allows moving from last entry if (source != numberOfBlocksP-1) ok = false; + if (source != numberOfBlocksR-1) ok = false; #ifdef USE_GPU if (cachedCapacity != currentCapacity) ok = false; #endif + if (currentCapacityR != currentCapacity) ok = false; + if (numberOfBlocksR != numberOfBlocks) ok = false; + if (currentCapacityG != currentCapacity) ok = false; + if (numberOfBlocksG != numberOfBlocks) ok = false; if (currentCapacityP != currentCapacity) ok = false; if (numberOfBlocksP != numberOfBlocks) ok = false; if (ok == false) { #if !defined(__CUDA_ARCH__) && !defined(__HIP_DEVICE_COMPILE__) std::stringstream ss; ss << "VBC ERROR: invalid source LID=" << source << " in copy, target=" << target << " #blocks=" << numberOfBlocks << " capacity=" << currentCapacity << std::endl; - ss << "or sizes are wrong, data->size()=" << block_data.size() << " parameters.size()=" << parameters.size() << std::endl; + ss << "or sizes are wrong, data->size()=" << block_data.size()<< " refined.size()=" << refined.size() << " ghost.size()=" << ghost.size() << " parameters.size()=" << parameters.size() << std::endl; std::cerr << ss.str(); sleep(1); exit(1); #else - printf("VBC error: invalid source LID=%u in copy, target=%u #blocks=%u capacity=%u \n or sizes are wrong, data->size()=%u parameters.size()=%u \n", - source,target,numberOfBlocks,currentCapacity, (vmesh::LocalID)block_data.size(),(vmesh::LocalID)parameters.size()); + printf("VBC error: invalid source LID=%u in copy, target=%u #blocks=%u capacity=%u \n or sizes are wrong, data->size()=%u refined.size()=%u ghost.size()=%u parameters.size()=%u \n", + source,target,numberOfBlocks,currentCapacity, (vmesh::LocalID)block_data.size(),(vmesh::LocalID)refined.size(),(vmesh::LocalID)ghost.size(),(vmesh::LocalID)parameters.size()); assert(0); #endif } @@ -289,12 +342,20 @@ namespace vmesh { for (unsigned int i=0; i= numberOfBlocks) { + exitInvalidLocalID(blockLID); + } + #else + if (blockLID >= numberOfBlocks) { + exitInvalidLocalID(blockLID,"getRefined"); + } + #endif + #endif + return refined.data() + blockLID*WID3; + } + + inline ARCH_HOSTDEV const uint8_t* VelocityBlockContainer::getRefined(const vmesh::LocalID blockLID) const { + #ifdef DEBUG_VBC + const vmesh::LocalID numberOfBlocks = refined.size()/WID3; + #if defined(USE_GPU) && (defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)) + if (blockLID >= numberOfBlocks) { + exitInvalidLocalID(blockLID); + } + #else + if (blockLID >= numberOfBlocks) { + exitInvalidLocalID(blockLID,"const getRefined const"); + } + #endif + #endif + return refined.data() + blockLID*WID3; + } + + inline ARCH_HOSTDEV uint8_t* VelocityBlockContainer::getGhost() { + return ghost.data(); + } + + inline ARCH_HOSTDEV const uint8_t* VelocityBlockContainer::getGhost() const { + return ghost.data(); + } + + inline ARCH_HOSTDEV uint8_t* VelocityBlockContainer::getGhost(const vmesh::LocalID blockLID) { + #ifdef DEBUG_VBC + const vmesh::LocalID numberOfBlocks = ghost.size(); + #if defined(USE_GPU) && (defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)) + if (blockLID >= numberOfBlocks) { + exitInvalidLocalID(blockLID); + } + #else + if (blockLID >= numberOfBlocks) { + exitInvalidLocalID(blockLID,"getGhost"); + } + #endif + #endif + return ghost.data() + blockLID; + } + + inline ARCH_HOSTDEV const uint8_t* VelocityBlockContainer::getGhost(const vmesh::LocalID blockLID) const { + #ifdef DEBUG_VBC + const vmesh::LocalID numberOfBlocks = ghost.size(); + #if defined(USE_GPU) && (defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)) + if (blockLID >= numberOfBlocks) { + exitInvalidLocalID(blockLID); + } + #else + if (blockLID >= numberOfBlocks) { + exitInvalidLocalID(blockLID,"const getGhost const"); + } + #endif + #endif + return ghost.data() + blockLID; + } + inline ARCH_HOSTDEV Real* VelocityBlockContainer::getParameters() { return parameters.data(); } @@ -429,6 +570,10 @@ namespace vmesh { } block_data.erase(block_data.begin() + WID3*(numberOfBlocks-1), block_data.begin() + WID3*(numberOfBlocks)); + refined.erase(refined.begin() + WID3*(numberOfBlocks-1), + refined.begin() + WID3*(numberOfBlocks)); + ghost.erase(ghost.begin() + (numberOfBlocks-1), + ghost.begin() + (numberOfBlocks)); parameters.erase(parameters.begin() + BlockParams::N_VELOCITY_BLOCK_PARAMS*(numberOfBlocks-1), parameters.begin() + BlockParams::N_VELOCITY_BLOCK_PARAMS*(numberOfBlocks)); #ifdef USE_GPU @@ -456,32 +601,40 @@ namespace vmesh { assert(0 && "ERROR! Attempting to grow block container on-device beyond capacity (::push_back)."); } block_data.device_resize((numberOfBlocks+1)*WID3); + ghost.device_resize(numberOfBlocks+1); + refined.device_resize((numberOfBlocks+1)*WID3); parameters.device_resize((numberOfBlocks+1)*BlockParams::N_VELOCITY_BLOCK_PARAMS); #elif defined(USE_GPU) setNewCapacity(numberOfBlocks+1,stream); block_data.resize((numberOfBlocks+1)*WID3,true,stream); + refined.resize((numberOfBlocks+1)*WID3,true,stream); + ghost.resize(numberOfBlocks+1,true,stream); parameters.resize((numberOfBlocks+1)*BlockParams::N_VELOCITY_BLOCK_PARAMS,true,stream); #else setNewCapacity(numberOfBlocks+1); block_data.resize((numberOfBlocks+1)*WID3,true); + refined.resize((numberOfBlocks+1)*WID3,true); + ghost.resize(numberOfBlocks+1,true); parameters.resize((numberOfBlocks+1)*BlockParams::N_VELOCITY_BLOCK_PARAMS,true); #endif #ifdef DEBUG_VBC const vmesh::LocalID currentCapacity = block_data.capacity()/WID3; + const vmesh::LocalID currentCapacityR = refined.capacity()/WID3; + const vmesh::LocalID currentCapacityG = ghost.capacity(); const vmesh::LocalID currentCapacityP = parameters.capacity()/BlockParams::N_VELOCITY_BLOCK_PARAMS; if (newIndex >= currentCapacity || newIndex >= currentCapacityP) { #if !defined(__CUDA_ARCH__) && !defined(__HIP_DEVICE_COMPILE__) std::stringstream ss; ss << "VBC ERROR in push_back, LID=" << newIndex << " for new block is out of bounds" << std::endl; - ss << "\t data->size()=" << block_data.size() << " parameters.size()=" << parameters.size() << std::endl; - ss << "\t data->capacity()=" << block_data.capacity() << " parameters.capacity()=" << parameters.capacity() << std::endl; + ss << "\t data->size()=" << block_data.size() << " refined.size()=" << refined.size() << " ghost.size()=" << ghost.size() << " parameters.size()=" << parameters.size() << std::endl; + ss << "\t data->capacity()=" << block_data.capacity() << " refined.capacity()=" << refined.capacity() << " ghost.capacity()=" << ghost.capacity() << " parameters.capacity()=" << parameters.capacity() << std::endl; std::cerr << ss.str(); sleep(1); exit(1); #else - printf("VBC ERROR in device push_back, LID=%u for new block is out of bounds\n data->size()=%u parameters.size()=%u\n", - newIndex,(vmesh::LocalID)block_data.size(),(vmesh::LocalID)parameters.size()); + printf("VBC ERROR in device push_back, LID=%u for new block is out of bounds\n data->size()=%u refined->size()=%u ghost->size()=%u parameters.size()=%u\n", + newIndex,(vmesh::LocalID)block_data.size(),(vmesh::LocalID)refined.size(),(vmesh::LocalID)ghost.size(),(vmesh::LocalID)parameters.size()); assert(0); #endif } @@ -513,32 +666,40 @@ namespace vmesh { assert(0 && "ERROR! Attempting to grow block container on-device beyond capacity (::push_back_and_zero)."); } block_data.device_resize((numberOfBlocks+1)*WID3, false); //construct=false don't construct or set to zero (performed below) + refined.device_resize((numberOfBlocks+1)*WID3, false); //construct=false don't construct or set to zero (performed below) + ghost.device_resize((numberOfBlocks+1), false); //construct=false don't construct or set to zero (performed below) parameters.device_resize((numberOfBlocks+1)*BlockParams::N_VELOCITY_BLOCK_PARAMS, false); //construct=false don't construct or set to zero (performed below) #elif defined(USE_GPU) setNewCapacity(numberOfBlocks+1,stream); block_data.resize((numberOfBlocks+1)*WID3,true,stream); + refined.resize((numberOfBlocks+1)*WID3,true,stream); + ghost.resize((numberOfBlocks+1),true,stream); parameters.resize((numberOfBlocks+1)*BlockParams::N_VELOCITY_BLOCK_PARAMS,true,stream); #else setNewCapacity(numberOfBlocks+1); block_data.resize((numberOfBlocks+1)*WID3); + refined.resize((numberOfBlocks+1)*WID3); + ghost.resize((numberOfBlocks+1)); parameters.resize((numberOfBlocks+1)*BlockParams::N_VELOCITY_BLOCK_PARAMS); #endif #ifdef DEBUG_VBC const vmesh::LocalID currentCapacity = block_data.capacity()/WID3; + const vmesh::LocalID currentCapacityR = refined.capacity()/WID3; + const vmesh::LocalID currentCapacityG = ghost.capacity(); const vmesh::LocalID currentCapacityP = parameters.capacity()/BlockParams::N_VELOCITY_BLOCK_PARAMS; - if (newIndex >= currentCapacity || newIndex >= currentCapacityP) { + if (newIndex >= currentCapacity || newIndex >= currentCapacityR || newIndex >= currentCapacityP) { #if !defined(__CUDA_ARCH__) && !defined(__HIP_DEVICE_COMPILE__) std::stringstream ss; ss << "VBC ERROR in push_back_and_zero, LID=" << newIndex << " for new block is out of bounds" << std::endl; - ss << "\t data->size()=" << block_data.size() << " parameters.size()=" << parameters.size() << std::endl; - ss << "\t data->capacity()=" << block_data.capacity() << " parameters.capacity()=" << parameters.capacity() << std::endl; + ss << "\t data->size()=" << block_data.size() << " refined.size()=" << refined.size() << " ghost.size()=" << ghost.size() << " parameters.size()=" << parameters.size() << std::endl; + ss << "\t data->capacity()=" << block_data.capacity() << " refined.capacity()=" << refined.capacity() << " ghost.capacity()=" << ghost.capacity() << " parameters.capacity()=" << parameters.capacity() << std::endl; std::cerr << ss.str(); sleep(1); exit(1); #else - printf("VBC ERROR in device push_back_and_zero, LID=%u for new block is out of bounds \n data->size()=%u parameters.size()=%u \n", - newIndex,(vmesh::LocalID)block_data.size(),(vmesh::LocalID)parameters.size()); + printf("VBC ERROR in device push_back_and_zero, LID=%u for new block is out of bounds \n data->size()=%u refined->size()=%u ghost->size()=%u parameters.size()=%u \n", + newIndex,(vmesh::LocalID)block_data.size(),(vmesh::LocalID)refined.size(),(vmesh::LocalID)ghost.size(),(vmesh::LocalID)parameters.size()); assert(0); #endif } @@ -547,6 +708,10 @@ namespace vmesh { for (size_t i=0; i> block_data_new(newCapacity*WID3); + std::vector refined_new(newCapacity*WID3); + std::vector ghost_new(newCapacity); std::vector> parameters_new(newCapacity*BlockParams::N_VELOCITY_BLOCK_PARAMS); block_data_new.resize(numberOfBlocks*WID3); + refined_new.resize(numberOfBlocks*WID3); + ghost_new.resize(numberOfBlocks); parameters_new.resize(numberOfBlocks*BlockParams::N_VELOCITY_BLOCK_PARAMS); for (size_t i=0; i= WID3) { + ok = false; + } + if (blockLID >= numberOfBlocks) { + ok = false; + } + if (blockLID*WID3+cell >= refined.size()) { + ok = false; + } + if (ok == false) { + std::stringstream ss; + ss << "VBC ERROR: out of bounds in getRefined, LID=" << blockLID << " cell=" << cell << " #blocks=" << numberOfBlocks << " refined->size()=" << refined.size() << std::endl; + std::cerr << ss.str(); + sleep(1); + exit(1); + } + return refined[blockLID*WID3+cell]; //May be bugged if it's not an integer + } + + inline const uint8_t& VelocityBlockContainer::getGhost(const vmesh::LocalID blockLID,const unsigned int cell) const { + const vmesh::LocalID numberOfBlocks = block_data.size()/WID3; + bool ok = true; + if (cell >= WID3) { + ok = false; + } + if (blockLID >= numberOfBlocks) { + ok = false; + } + if (blockLID+cell >= ghost.size()) { + ok = false; + } + if (ok == false) { + std::stringstream ss; + ss << "VBC ERROR: out of bounds in getGhost, LID=" << blockLID << " cell=" << cell << " #blocks=" << numberOfBlocks << " ghost->size()=" << ghost.size() << std::endl; + std::cerr << ss.str(); + sleep(1); + exit(1); + } + return ghost[blockLID+cell]; //May be bugged if it's not an integer + } + inline const Real& VelocityBlockContainer::getParameters(const vmesh::LocalID blockLID,const unsigned int cell) const { const vmesh::LocalID numberOfBlocks = block_data.size()/WID3; bool ok = true; @@ -905,6 +1160,53 @@ namespace vmesh { block_data[blockLID*WID3+cell] = value; } + + inline void VelocityBlockContainer::setRefined(const vmesh::LocalID blockLID,const unsigned int cell,const uint8_t value) { + const vmesh::LocalID numberOfBlocks = block_data.size()/WID3; + bool ok = true; + if (cell >= WID3) { + ok = false; + } + if (blockLID >= numberOfBlocks) { + ok = false; + } + if (blockLID*WID3+cell >= refined.size()) { + ok = false; + } + if (ok == false) { + std::stringstream ss; + ss << "VBC ERROR: out of bounds in setRefined, LID=" << blockLID << " cell=" << cell << " #blocks=" << numberOfBlocks << " refined->size()=" << refined.size() << std::endl; + std::cerr << ss.str(); + sleep(1); + exit(1); + } + + refined[blockLID*WID3+cell] = value; + } + + inline void VelocityBlockContainer::setGhost(const vmesh::LocalID blockLID,const unsigned int cell,const uint8_t value) { + const vmesh::LocalID numberOfBlocks = block_data.size()/WID3; + bool ok = true; + if (cell >= WID3) { + ok = false; + } + if (blockLID >= numberOfBlocks) { + ok = false; + } + if (blockLID+cell >= ghost.size()) { + ok = false; + } + if (ok == false) { + std::stringstream ss; + ss << "VBC ERROR: out of bounds in setGhost, LID=" << blockLID << " cell=" << cell << " #blocks=" << numberOfBlocks << " ghost->size()=" << ghost.size() << std::endl; + std::cerr << ss.str(); + sleep(1); + exit(1); + } + + ghost[blockLID+cell] = value; + } + #endif //debug VBC } // namespace block_cont diff --git a/vlasiator.cpp b/vlasiator.cpp index acbb9ad59f..13ed6b5397 100644 --- a/vlasiator.cpp +++ b/vlasiator.cpp @@ -670,7 +670,7 @@ int simulate(int argn,char* args[]) { // Run Vlasov solver once with zero dt to initialize // per-cell dt limits. Also compute initial _R and _V moments at restart. calculateSpatialTranslation(mpiGrid,0.0); - calculateAcceleration(mpiGrid,0.0); + calculateAcceleration(mpiGrid,0.0,true); sysBoundaryContainer.setupL2OutflowAtRestart(mpiGrid); @@ -741,10 +741,10 @@ int simulate(int argn,char* args[]) { //the distribution function is already propagated forward in time by dt/2 phiprof::Timer propagateHalfTimer {"propagate-velocity-space-dt/2"}; if (P::propagateVlasovAcceleration) { - calculateAcceleration(mpiGrid, 0.5*P::dt); + calculateAcceleration(mpiGrid, 0.5*P::dt,true); } else { //zero step to set up moments _v - calculateAcceleration(mpiGrid, 0.0); + calculateAcceleration(mpiGrid, 0.0,true); } propagateHalfTimer.stop(); @@ -1121,7 +1121,7 @@ int simulate(int argn,char* args[]) { // Calculate new dt limits since we might break CFL when refining phiprof::Timer computeDtimer {"compute-dt-amr"}; calculateSpatialTranslation(mpiGrid,0.0); - calculateAcceleration(mpiGrid,0.0); + calculateAcceleration(mpiGrid,0.0,true); } // This now uses the block-based count just copied between the two refinement calls above. balanceLoad(mpiGrid, sysBoundaryContainer, technical.view(), fsgrid); @@ -1158,11 +1158,11 @@ int simulate(int argn,char* args[]) { //propagate velocity space back to real-time if( P::propagateVlasovAcceleration ) { // Back half dt to real time, forward by new half dt - calculateAcceleration(mpiGrid,-0.5*P::dt + 0.5*newDt); + calculateAcceleration(mpiGrid,-0.5*P::dt + 0.5*newDt,false); } else { //zero step to set up moments _v - calculateAcceleration(mpiGrid, 0.0); + calculateAcceleration(mpiGrid, 0.0,false); } P::dt=newDt; @@ -1331,11 +1331,19 @@ int simulate(int argn,char* args[]) { phiprof::Timer vspaceTimer {"Velocity-space"}; if ( P::propagateVlasovAcceleration ) { - calculateAcceleration(mpiGrid,P::dt); + if(P::activateVamr && (P::tstep % P::vAMRrefineStep) == 0){ + calculateAcceleration(mpiGrid,P::dt, true); + }else{ + calculateAcceleration(mpiGrid,P::dt, false); + } addTimedBarrier("barrier-after-ad just-blocks"); } else { //zero step to set up moments _v - calculateAcceleration(mpiGrid, 0.0); + if(P::activateVamr && (P::tstep % P::vAMRrefineStep) == 0){ + calculateAcceleration(mpiGrid, 0.0,true); + }else{ + calculateAcceleration(mpiGrid, 0.0,false); + } } vspaceTimer.stop(computedCells, "Cells"); addTimedBarrier("barrier-after-acceleration"); diff --git a/vlasovsolver/arch_moments.cpp b/vlasovsolver/arch_moments.cpp index 296bc03b69..53de9639ab 100644 --- a/vlasovsolver/arch_moments.cpp +++ b/vlasovsolver/arch_moments.cpp @@ -85,21 +85,59 @@ void calculateCellMoments(spatial_cell::SpatialCell* cell, continue; } + vmesh::MeshParameters& vMeshprint=vmesh::getMeshWrapper()->velocityMeshesCreation->at(popID); + species::Species& species=getObjectWrapper().particleSpecies[popID]; + const Real mass = getObjectWrapper().particleSpecies[popID].mass; const Real charge = getObjectWrapper().particleSpecies[popID].charge; - // Temporary array for storing moments + // Temporary array where the moments for this species are accumulated Real array[nMom1] = {0}; + + int Ref=0; + int MaxRef=0; + if(P::activateVamr) { + Ref=getObjectWrapper().particleSpecies[popID].RefinementLevel; + MaxRef=getObjectWrapper().particleSpecies[popID].MaxRefinementLevel; + if(Ref < MaxRef){ + // Update the Refined parameter to know the cells that need to be integrated + changeRefined(cell,popID); + } + // Calculate species' contribution to first velocity moments with vAMR + blockVelocityFirstMomentsVamr(blockContainer, + array, + nBlocks); + }else { + // Calculate species' contribution to first velocity moments + blockVelocityFirstMoments(blockContainer, + array, + nBlocks); + } + - // Calculate species' contribution to first velocity moments - blockVelocityFirstMoments(blockContainer, - array, - nBlocks); pop.RHO = array[0]; pop.V[0] = divideIfNonZero(array[1], array[0]); pop.V[1] = divideIfNonZero(array[2], array[0]); pop.V[2] = divideIfNonZero(array[3], array[0]); + if(P::activateVamr && Ref==MaxRef && MaxRef>0 ){ + //Sum of all the partial integrals and sharing of the final result between the different vAMR grids of the same species + for (uint popID2=(popID-MaxRef); popID2get_population(popID2); + pop.RHO += pop2.RHO; + pop.V[0] += pop2.V[0]; + pop.V[1] += pop2.V[1]; + pop.V[2] += pop2.V[2]; + }; + for (uint popID2=(popID-MaxRef); popID2get_population(popID2); + pop2.RHO = pop.RHO; + pop2.V[0] = pop.V[0]; + pop2.V[1] = pop.V[1]; + pop2.V[2] = pop.V[2]; + }; + }; + if (!computePopulationMomentsOnly) { // Store species' contribution to bulk velocity moments cell->parameters[CellParams::RHOM ] += array[0]*mass; @@ -140,13 +178,27 @@ void calculateCellMoments(spatial_cell::SpatialCell* cell, // Temporary array for storing moments Real array[nMom2] = {0}; - // Calculate species' contribution to second velocity moments - blockVelocitySecondMoments(blockContainer, - cell->parameters[CellParams::VX], - cell->parameters[CellParams::VY], - cell->parameters[CellParams::VZ], - array, - nBlocks); + int Ref=0; + int MaxRef=0; + if(P::activateVamr) { + int Ref=getObjectWrapper().particleSpecies[popID].RefinementLevel; + int MaxRef=getObjectWrapper().particleSpecies[popID].MaxRefinementLevel; + // Calculate species' contribution to second velocity moments with vAMR + blockVelocitySecondMomentsVamr(blockContainer, + cell->parameters[CellParams::VX], + cell->parameters[CellParams::VY], + cell->parameters[CellParams::VZ], + array, + nBlocks); + }else{ + // Calculate species' contribution to second velocity moments + blockVelocitySecondMoments(blockContainer, + cell->parameters[CellParams::VX], + cell->parameters[CellParams::VY], + cell->parameters[CellParams::VZ], + array, + nBlocks); + } // Store species' contribution to bulk velocity moments Population &pop = cell->get_population(popID); for (size_t i=0; iparameters[CellParams::P_13] += pop.P[4]; cell->parameters[CellParams::P_12] += pop.P[5]; } + + if(P::activateVamr && Ref==MaxRef && MaxRef>0 ){ + //Sum of all the partial integrals and sharing of the final result between the different vAMR grids of the same species + for (uint popID2=(popID-MaxRef); popID2get_population(popID2); + for (size_t i=0; iget_population(popID2); + for (size_t i=0; i0 ){ + //Sum of all the partial integrals and sharing of the final result between the different vAMR grids of the same species + for (uint popID2=(popID-MaxRef); popID2get_population(popID2); + pop.RHO_R += pop2.RHO_R; + pop.V_R[0] += pop2.V_R[0]; + pop.V_R[1] += pop2.V_R[1]; + pop.V_R[2] += pop2.V_R[2]; + }; + for (uint popID2=(popID-MaxRef); popID2get_population(popID2); + pop2.RHO_R = pop.RHO_R; + pop2.V_R[0] = pop.V_R[0]; + pop2.V_R[1] = pop.V_R[1]; + pop2.V_R[2] = pop.V_R[2]; + }; + }; + cell->parameters[CellParams::RHOM_R ] += array[0]*mass; cell->parameters[CellParams::VX_R] += array[1]*mass; cell->parameters[CellParams::VY_R] += array[2]*mass; @@ -305,13 +408,28 @@ void calculateMoments_R( // Temporary array where species' contribution to 2nd moments is accumulated Real array[nMom2] = {0}; - // Calculate species' contribution to second velocity moments - blockVelocitySecondMoments(blockContainer, - cell->parameters[CellParams::VX_R], - cell->parameters[CellParams::VY_R], - cell->parameters[CellParams::VZ_R], - array, - nBlocks); + int Ref=0; + int MaxRef=0; + if(P::activateVamr) { + int Ref=getObjectWrapper().particleSpecies[popID].RefinementLevel; + int MaxRef=getObjectWrapper().particleSpecies[popID].MaxRefinementLevel; + // Calculate species' contribution to second velocity moments with Vam + blockVelocitySecondMomentsVamr(blockContainer, + cell->parameters[CellParams::VX], + cell->parameters[CellParams::VY], + cell->parameters[CellParams::VZ], + array, + nBlocks); + }else{ + // Calculate species' contribution to second velocity moments + blockVelocitySecondMoments(blockContainer, + cell->parameters[CellParams::VX], + cell->parameters[CellParams::VY], + cell->parameters[CellParams::VZ], + array, + nBlocks); + } + // Store species' contribution to 2nd bulk velocity moments Population &pop = cell->get_population(popID); for (size_t i = 0; i < nMom2; ++i) { @@ -324,6 +442,23 @@ void calculateMoments_R( cell->parameters[CellParams::P_23_R] += pop.P_R[3]; cell->parameters[CellParams::P_13_R] += pop.P_R[4]; cell->parameters[CellParams::P_12_R] += pop.P_R[5]; + + if(P::activateVamr && Ref==MaxRef && MaxRef>0 ){ + //Sum of all the partial integrals and sharing of the final result between the different vAMR grids of the same species + for (uint popID2=(popID-MaxRef); popID2get_population(popID2); + for (size_t i=0; iget_population(popID2); + for (size_t i=0; i0 ){ + //Sum of all the partial integrals and sharing of the final result between the different vAMR grids of the same species + for (uint popID2=(popID-MaxRef); popID2get_population(popID2); + pop.RHO_V += pop2.RHO_V; + pop.V_V[0] += pop2.V_V[0]; + pop.V_V[1] += pop2.V_V[1]; + pop.V_V[2] += pop2.V_V[2]; + }; + for (uint popID2=(popID-MaxRef); popID2get_population(popID2); + pop2.RHO_V = pop.RHO_V; + pop2.V_V[0] = pop.V_V[0]; + pop2.V_V[1] = pop.V_V[1]; + pop2.V_V[2] = pop.V_V[2]; + }; + }; cell->parameters[CellParams::RHOM_V ] += array[0]*mass; cell->parameters[CellParams::VX_V] += array[1]*mass; @@ -475,13 +644,28 @@ void calculateMoments_V( // Temporary array where moments are stored Real array[nMom2] = {0}; - // Calculate species' contribution to second velocity moments - blockVelocitySecondMoments(blockContainer, - cell->parameters[CellParams::VX_V], - cell->parameters[CellParams::VY_V], - cell->parameters[CellParams::VZ_V], - array, - nBlocks); + int Ref=0; + int MaxRef=0; + if(P::activateVamr) { + int Ref=getObjectWrapper().particleSpecies[popID].RefinementLevel; + int MaxRef=getObjectWrapper().particleSpecies[popID].MaxRefinementLevel; + // Calculate species' contribution to second velocity moments with vAMR + blockVelocitySecondMomentsVamr(blockContainer, + cell->parameters[CellParams::VX], + cell->parameters[CellParams::VY], + cell->parameters[CellParams::VZ], + array, + nBlocks); + }else{ + // Calculate species' contribution to second velocity moments + blockVelocitySecondMoments(blockContainer, + cell->parameters[CellParams::VX], + cell->parameters[CellParams::VY], + cell->parameters[CellParams::VZ], + array, + nBlocks); + } + // Store species' contribution to 2nd bulk velocity moments Population &pop = cell->get_population(popID); for (size_t i = 0; i < nMom2; ++i) { @@ -494,6 +678,1143 @@ void calculateMoments_V( cell->parameters[CellParams::P_23_V] += pop.P_V[3]; cell->parameters[CellParams::P_13_V] += pop.P_V[4]; cell->parameters[CellParams::P_12_V] += pop.P_V[5]; + + if(P::activateVamr && Ref==MaxRef && MaxRef>0 ){ + //Sum of all the partial integrals and sharing of the final result between the different vAMR grids of the same species + for (uint popID2=(popID-MaxRef); popID2get_population(popID2); + for (size_t i=0; iget_population(popID2); + for (size_t i=0; i& mpiGrid, + const std::vector& cells, + const uint popID) { + + // Loop over all particle species +#pragma omp parallel for schedule(dynamic,1) + for (size_t c=0; cgetVelocityBlockMinValue(popID); + +#ifdef USE_GPU + vmesh::VelocityMesh* vmesh = cell->dev_get_velocity_mesh(popID); + vmesh::VelocityMesh* vmeshraf = cell->get_velocity_mesh(popID+1); + vmesh::VelocityBlockContainer* blockContainer = cell->dev_get_velocity_blocks(popID); + vmesh::VelocityBlockContainer* blockContainerraf = cell->dev_get_velocity_blocks(popID+1); +#else + vmesh::VelocityMesh* vmesh = cell->get_velocity_mesh(popID); + vmesh::VelocityMesh* vmeshraf = cell->get_velocity_mesh(popID+1); + vmesh::VelocityBlockContainer* blockContainer = cell->get_velocity_blocks(popID); + vmesh::VelocityBlockContainer* blockContainerraf = cell->get_velocity_blocks(popID+1); +#endif + const uint nBlocks = cell->get_velocity_mesh(popID)->size(); + Population &pop = cell->get_population(popID); + if (nBlocks == 0) { + continue; + } + + Realf *data = blockContainer->getData(); + Realf *dataraf = blockContainerraf->getData(); + uint8_t *ghost = blockContainerraf->getGhost(); + + for (vmesh::LocalID localID=0; localIDsize(); ++localID) { + const vmesh::GlobalID globalID = vmesh->getGlobalID(localID); + + vmesh::LocalID Indices[3]; + vmesh->getIndices(globalID, Indices[0], Indices[1], Indices[2]); + + for (int i=0; i<2; ++i) { + for (int j=0; j<2; ++j) { + for (int k=0; k<2; ++k) { + //Indices of the refined blocks + //Each coarse block contains 8 refined blocks + vmesh::LocalID Indicesraf[3]; + Indicesraf[0] = 2*Indices[0]+i ; + Indicesraf[1] = 2*Indices[1]+j ; + Indicesraf[2] = 2*Indices[2]+k ; + + const vmesh::GlobalID globalIDraf=vmeshraf->getGlobalID(Indicesraf); + //Check if the refined cell exists + if (globalIDraf == vmeshraf->invalidGlobalID()) { + continue; + }else{ + const vmesh::LocalID localIDraf=vmeshraf->getLocalID(globalIDraf); + if (localIDraf == vmeshraf->invalidLocalID()) { + continue; + }else{ + for (int i2=0; i2<2; ++i2) { + for (int j2=0; j2<2; ++j2) { + for (int k2=0; k2<2; ++k2) { + //Loop over the 8 coarsed cells located on the WID3 refined cells (the refined block) + Realf summ=0; + // if( data[localID*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]>cell->getVelocityBlockMinValue(popID+1)){ //old criteria used to improve communication + Realf datasave= data[localID*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]; + data[localID*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]=0; + for (int i3=0; i3<2; ++i3) { + for (int j3=0; j3<2; ++j3) { + for (int k3=0; k3<2; ++k3) { + //Loop over the 8 refined cells contained in the 1 coarsed cell + if(dataraf[localIDraf*WID3+cellIndex(2*i2+i3,2*j2+j3,2*k2+k3)]>minValue){ + // if(ghost[localIDraf]==1 && dataraf[localIDraf*WID3+cellIndex(2*i2+i3,2*j2+j3,2*k2+k3)]>minValue){ //Another criteria + data[localID*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]+= dataraf[localIDraf*WID3+cellIndex(2*i2+i3,2*j2+j3,2*k2+k3)]/8.0; + summ+=1.0; + }else{ + dataraf[localIDraf*WID3+cellIndex(2*i2+i3,2*j2+j3,2*k2+k3)]=datasave; + } + } + } + } + if (summ!=8.0){ + data[localID*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]=datasave; + } + /*}else{ + for (int i3=0; i3<2; ++i3) { + for (int j3=0; j3<2; ++j3) { + for (int k3=0; k3<2; ++k3) { + dataraf[localIDraf*WID3+cellIndex(2*i2+i3,2*j2+j3,2*k2+k3)]=data[localID*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]; + } + } + } + }*/ //Removed because the boundaries were looking too coarse + + }; + }; + };//loop over refined cells + }; + }; + };//loop over refined blocks + }; + }; + }; + } +} + +// Update the Refined parameter to know the cells that need to be integrated + void changeRefined(spatial_cell::SpatialCell* cell, + const uint popID){ + + vmesh::VelocityMesh* vmesh = cell->get_velocity_mesh(popID); + vmesh::VelocityMesh* vmeshraf = cell->get_velocity_mesh(popID+1); + + uint8_t *refined =cell->get_velocity_blocks(popID)->getRefined(); +//Parallelisation is just in case; this function may already be part of a parallelised region +#pragma omp parallel for schedule(dynamic,1) + for (vmesh::LocalID localID=0; localIDsize(); ++localID) { + const vmesh::GlobalID globalID = vmesh->getGlobalID(localID); + + vmesh::LocalID Indices[3]; + vmesh->getIndices(globalID, Indices[0], Indices[1], Indices[2]); + + for (int i=0; i<2; ++i) { + for (int j=0; j<2; ++j) { + for (int k=0; k<2; ++k) { + //Indices of the refined blocks + //Each coarse block contains 8 refined blocks + vmesh::LocalID Indicesraf[3]; + Indicesraf[0] = 2*Indices[0]+i ; + Indicesraf[1] = 2*Indices[1]+j ; + Indicesraf[2] = 2*Indices[2]+k ; + + const vmesh::GlobalID globalIDraf=vmeshraf->getGlobalID(Indicesraf); + //We will check if the refined block exists + if (globalIDraf == vmeshraf->invalidGlobalID()) { + //each refined block represent (WID/2)^3 coarse cells (8 in our case) + refined[localID*WID3+cellIndex(2*i,2*j,2*k)]=false; + refined[localID*WID3+cellIndex(2*i+1,2*j,2*k)]=false; + refined[localID*WID3+cellIndex(2*i,2*j+1,2*k)]=false; + refined[localID*WID3+cellIndex(2*i,2*j,2*k+1)]=false; + refined[localID*WID3+cellIndex(2*i+1,2*j+1,2*k)]=false; + refined[localID*WID3+cellIndex(2*i,2*j+1,2*k+1)]=false; + refined[localID*WID3+cellIndex(2*i+1,2*j,2*k+1)]=false; + refined[localID*WID3+cellIndex(2*i+1,2*j+1,2*k+1)]=false; + }else{ + const vmesh::LocalID localIDraf=vmeshraf->getLocalID(globalIDraf); + if (localIDraf == vmeshraf->invalidLocalID()) { + refined[localID*WID3+cellIndex(2*i,2*j,2*k)]=false; + refined[localID*WID3+cellIndex(2*i+1,2*j,2*k)]=false; + refined[localID*WID3+cellIndex(2*i,2*j+1,2*k)]=false; + refined[localID*WID3+cellIndex(2*i,2*j,2*k+1)]=false; + refined[localID*WID3+cellIndex(2*i+1,2*j+1,2*k)]=false; + refined[localID*WID3+cellIndex(2*i,2*j+1,2*k+1)]=false; + refined[localID*WID3+cellIndex(2*i+1,2*j,2*k+1)]=false; + refined[localID*WID3+cellIndex(2*i+1,2*j+1,2*k+1)]=false; + }else{ + //The refine block cell + refined[localID*WID3+cellIndex(2*i,2*j,2*k)]=true; + refined[localID*WID3+cellIndex(2*i+1,2*j,2*k)]=true; + refined[localID*WID3+cellIndex(2*i,2*j+1,2*k)]=true; + refined[localID*WID3+cellIndex(2*i,2*j,2*k+1)]=true; + refined[localID*WID3+cellIndex(2*i+1,2*j+1,2*k)]=true; + refined[localID*WID3+cellIndex(2*i,2*j+1,2*k+1)]=true; + refined[localID*WID3+cellIndex(2*i+1,2*j,2*k+1)]=true; + refined[localID*WID3+cellIndex(2*i+1,2*j+1,2*k+1)]=true; + }; + }; + }; + }; + }; + }; + } + +//Every vamr_refinedStep we check all the velocity cells with the vAMR criterion with the 1st order +void RefinedOrder1(dccrg::Dccrg& mpiGrid, + const std::vector& cells){ + +#pragma omp parallel for schedule(dynamic,1) +for (size_t c=0; c ListBlockExist[getObjectWrapper().particleSpecies.size()]; + + for (int popID=(getObjectWrapper().particleSpecies.size()-2); popID>-1; --popID) { + + if(getObjectWrapper().particleSpecies[popID].MaxRefinementLevel>0 && getObjectWrapper().particleSpecies[popID].RefinementLevelget_velocity_mesh(popID); + vmesh::VelocityMesh* vmeshraf = cell->get_velocity_mesh(popID+1); + Realf *data = cell->get_velocity_blocks(popID)->getData(); + + for (vmesh::LocalID localID=0; localIDsize(); ++localID) { + vmesh::GlobalID globalID = vmesh->getGlobalID(localID); + vmesh::LocalID Indices[3]; + vmesh->getIndices(globalID, Indices[0], Indices[1], Indices[2]); + + if (getObjectWrapper().particleSpecies[popID].RefinementLevel==0) { + ListBlockExist[popID].insert(globalID); + } + + for (int i=0; i<2; ++i) { + for (int j=0; j<2; ++j) { + for (int k=0; k<2; ++k) { + //Loop over the future refined blocks R+1 + Realf Datagros = 0; + + for (int i2=0; i2<2; ++i2) { + for (int j2=0; j2<2; ++j2) { + for (int k2=0; k2<2; ++k2) { + //Sum over the 8 cells of level R in order to reproduce the coarse cell R-1 + Datagros += data[localID*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]; + } + } + } + Datagros/=8.0; + + Realf D = abs( data[localID*WID3+cellIndex(1+i,1+j,1+k)] - Datagros ); + // The idea is to always compare the 8 central cells of level R with the reproduce R-1 cells + + vmesh::LocalID Indicesraf[3]; + Indicesraf[0] = 2*Indices[0]+i; + Indicesraf[1] = 2*Indices[1]+j; + Indicesraf[2] = 2*Indices[2]+k; + + Realf Dcomp = cell->getVelocityBlockMinValue(popID); + + if(getObjectWrapper().particleSpecies[popID].CriteriaMethod==1){ //epsilon + Dcomp= getObjectWrapper().particleSpecies[popID].CriteriaValue; + }else if(getObjectWrapper().particleSpecies[popID].CriteriaMethod==2){ //epsilon*2^-R + Dcomp= getObjectWrapper().particleSpecies[popID].CriteriaValue/(1u << getObjectWrapper().particleSpecies[popID].RefinementLevel); + } + if (D > Dcomp){ + // We should create a new block for R+1 + int addWidthV = 1; //getObjectWrapper().particleSpecies[popID+1].sparseBlockAddWidthV; + for (int offset_vx=-addWidthV;offset_vx<=addWidthV;offset_vx++) { + for (int offset_vy=-addWidthV;offset_vy<=addWidthV;offset_vy++) { + for (int offset_vz=-addWidthV;offset_vz<=addWidthV;offset_vz++) { + const vmesh::GlobalID globalIDraf = vmeshraf->getGlobalID(Indicesraf[0]+offset_vx,Indicesraf[1]+offset_vy,Indicesraf[2]+offset_vz); + if (globalIDraf== vmeshraf->invalidGlobalID()) { + // std::cout<< " GlobalID bug not normal" << "Indices[0]+offset_vx" << Indices[0]+offset_vx << "Indices[1]+offset_vy " << Indices[1]+offset_vy << "Indices[1]+offset_vy" << Indices[2]+offset_vz <-1; --popID) { + + if(getObjectWrapper().particleSpecies[popID].MaxRefinementLevel>0 && getObjectWrapper().particleSpecies[popID].RefinementLevelget_velocity_mesh(popID); + vmesh::VelocityMesh* vmeshraf = cell->get_velocity_mesh(popID+1); + + for (vmesh::GlobalID globalIDraf : ListBlockExist[popID+1]) { + + vmesh::LocalID Indicesraf[3]; + vmeshraf->getIndices(globalIDraf, Indicesraf[0], Indicesraf[1], Indicesraf[2]); + + vmesh::LocalID Indices[3]; + Indices[0] = Indicesraf[0]/2; + Indices[1] = Indicesraf[1]/2; + Indices[2] = Indicesraf[2]/2; + + int addWidthV = 0; //getObjectWrapper().particleSpecies[popID].sparseBlockAddWidthV; //Replace species.sparseBlockAddWidthV + for (int offset_vx=-addWidthV;offset_vx<=addWidthV;offset_vx++) { + for (int offset_vy=-addWidthV;offset_vy<=addWidthV;offset_vy++) { + for (int offset_vz=-addWidthV;offset_vz<=addWidthV;offset_vz++) { + const vmesh::GlobalID globalID = vmesh->getGlobalID(Indices[0]+offset_vx,Indices[1]+offset_vy,Indices[2]+offset_vz); + ListBlockExist[popID].insert(globalID); + } + } + } + } + + } + } + + //create the cells if needed + for (uint popID=0; popID0 ){ + + vmesh::VelocityMesh* vmesh = cell->get_velocity_mesh(popID); + vmesh::LocalID Localsize= vmesh->size(); + + for (vmesh::LocalID localID=0; localIDget_velocity_mesh(popID); + vmesh::GlobalID globalID = vmesh->getGlobalID(localID); + if (ListBlockExist[popID].find(globalID) == ListBlockExist[popID].end()) { + cell->remove_velocity_block(globalID,popID); + Localsize-=1; + localID-=1; + } + } + + + if(getObjectWrapper().particleSpecies[popID].RefinementLevel!=0){ + + for (vmesh::GlobalID globalID : ListBlockExist[popID]) { + + if(cell->add_velocity_block(globalID,popID)){ //True if it's a new block + //need to adapt the creation to the level of refinement needed + vmesh::VelocityMesh* vmesh = cell->get_velocity_mesh(popID); + vmesh::VelocityMesh* vmeshgros= cell->get_velocity_mesh(popID-1); + Realf *datagros = cell->get_velocity_blocks(popID-1)->getData(); + Realf *data = cell->get_velocity_blocks(popID)->getData(); + + vmesh::LocalID Indices[3]; + vmesh->getIndices(globalID, Indices[0], Indices[1], Indices[2]); + // i, j and k indicate the separation of the coarse block into eight refined blocks + int i = Indices[0]%2; + int j = Indices[1]%2; + int k = Indices[2]%2; + + vmesh::LocalID Indicesgros[3]; + Indicesgros[0] = (Indices[0]-i)/2; + Indicesgros[1] = (Indices[1]-j)/2; + Indicesgros[2] = (Indices[2]-k)/2; + + vmesh::GlobalID globalIDgros=vmeshgros->getGlobalID(Indicesgros); + vmesh::LocalID localIDgros=vmeshgros->getLocalID(globalIDgros); + vmesh::LocalID localIDcreated=vmesh->getLocalID(globalID); + uint8_t *ghost = cell->get_velocity_blocks(popID)->getGhost(localIDcreated); + ghost[0]=1; + + for (int i2=0; i2<2; ++i2) { + for (int j2=0; j2<2; ++j2) { + for (int k2=0; k2<2; ++k2) { + //Loop over the 8 coarsed cells located on the WID3 refined cells (the refined block) + data[localIDcreated*WID3+cellIndex(2*i2,2*j2,2*k2)]=datagros[localIDgros*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]; + data[localIDcreated*WID3+cellIndex(2*i2+1,2*j2,2*k2)]=datagros[localIDgros*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]; + data[localIDcreated*WID3+cellIndex(2*i2,2*j2+1,2*k2)]=datagros[localIDgros*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]; + data[localIDcreated*WID3+cellIndex(2*i2,2*j2,2*k2+1)]=datagros[localIDgros*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]; + data[localIDcreated*WID3+cellIndex(2*i2+1,2*j2+1,2*k2)]=datagros[localIDgros*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]; + data[localIDcreated*WID3+cellIndex(2*i2+1,2*j2,2*k2+1)]=datagros[localIDgros*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]; + data[localIDcreated*WID3+cellIndex(2*i2,2*j2+1,2*k2+1)]=datagros[localIDgros*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]; + data[localIDcreated*WID3+cellIndex(2*i2+1,2*j2+1,2*k2+1)]=datagros[localIDgros*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]; + } + } + } + + }else{ + vmesh::VelocityMesh* vmesh = cell->get_velocity_mesh(popID); + vmesh::LocalID localIDcreated=vmesh->getLocalID(globalID); + uint8_t *ghost = cell->get_velocity_blocks(popID)->getGhost(localIDcreated); + ghost[0]=1; + } + } + }else{ + for (vmesh::GlobalID globalID : ListBlockExist[popID]) { + cell->add_velocity_block(globalID,popID); + //ghost not important for this one + vmesh::VelocityMesh* vmesh = cell->get_velocity_mesh(popID); + vmesh::LocalID localIDcreated=vmesh->getLocalID(globalID); + uint8_t *ghost = cell->get_velocity_blocks(popID)->getGhost(localIDcreated); + ghost[0]=1; + } + } + + } + } +} +} + +//Every vamr_refinedStep we check all the velocity cells with the vAMR criterion with the 3rd order +void RefinedOrder3(dccrg::Dccrg& mpiGrid, + const std::vector& cells){ + + Realf M[3][3][3]; + Realf A[3] = {1.0/8.0, 1.0, -1.0/8.0}; + + for (int i0 = 0; i0 < 3; i0++) { + for (int j0 = 0; j0 < 3; j0++) { + for (int k0 = 0; k0 < 3; k0++) { + M[i0][j0][k0] = A[i0] * A[j0] * A[k0]; + } + } + } + +#pragma omp parallel for schedule(dynamic,1) +for (size_t c=0; c ListBlockExist[getObjectWrapper().particleSpecies.size()]; + + for (int popID=(getObjectWrapper().particleSpecies.size()-2); popID>-1; --popID) { + + if(getObjectWrapper().particleSpecies[popID].MaxRefinementLevel>0 && getObjectWrapper().particleSpecies[popID].RefinementLevelget_velocity_mesh(popID); + vmesh::VelocityMesh* vmeshraf = cell->get_velocity_mesh(popID+1); + Realf *data = cell->get_velocity_blocks(popID)->getData(); + + for (vmesh::LocalID localID=0; localIDsize(); ++localID) { + vmesh::GlobalID globalID = vmesh->getGlobalID(localID); + vmesh::LocalID Indices[3]; + vmesh->getIndices(globalID, Indices[0], Indices[1], Indices[2]); + + if (getObjectWrapper().particleSpecies[popID].RefinementLevel==0) { + ListBlockExist[popID].insert(globalID); + } + + bool Verif[8]; + for (int fill=0; fill<8; ++fill){ + Verif[fill]=true; + } + //We will need 1 ghost cell of level R-1 in every direction, meaning that we need a 4x4x4 grid of level R-1 cells for our block of level R + Realf Datagros[4][4][4]; + + for (int neighbour_vx=0; neighbour_vx<4; ++neighbour_vx) { + int n_vx=(neighbour_vx+1)/2 -1; //gives the shift in block + int i=abs((neighbour_vx+1)%2); + for (int neighbour_vy=0; neighbour_vy<4; ++neighbour_vy) { + int n_vy=(neighbour_vy+1)/2 -1; + int j=abs((neighbour_vy+1)%2); + for (int neighbour_vz=0; neighbour_vz<4; ++neighbour_vz) { + int n_vz=(neighbour_vz+1)/2 -1; + int k=abs((neighbour_vz+1)%2); + + vmesh::LocalID localIDneigh = vmesh->invalidLocalID(); + + if(n_vx==0 && n_vy==0 && n_vz==0){ + //The cell of level R-1 will be defined on the actual block of level R + localIDneigh = localID; + }else{ + //The cell of level R-1 will be defined on a neighbouring block of level R + const vmesh::GlobalID globalIDneigh = vmesh->getGlobalID(Indices[0]+n_vx,Indices[1]+n_vy,Indices[2]+n_vz); + if (globalIDneigh != vmesh->invalidGlobalID()) { + localIDneigh=vmesh->getLocalID(globalIDneigh); + } + } + + if (localIDneigh != vmesh->invalidLocalID()) { + //The block exists, meaning that we can compute the cell value for level R-1 + Realf Datagrossum=0; + for (int i2=0; i2<2; ++i2) { + for (int j2=0; j2<2; ++j2) { + for (int k2=0; k2<2; ++k2) { + Datagrossum += data[localIDneigh*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]; + } + } + } + //Computation of the grid R-1 + Datagros[neighbour_vx][neighbour_vy][neighbour_vz] = Datagrossum/8.0; + }else{ + for(int iloop = 0; iloop< (2-abs(n_vx)); ++iloop){ + for(int jloop = 0; jloop< (2-abs(n_vy)); ++jloop){ + for(int kloop = 0; kloop< (2-abs(n_vz)); ++kloop){ + //Loop over all the future refined blocks that will be impacted by the non-existing block + //It is technically impossible for the block in which we have the local ID to return invalidLocalID(), but this has been taken into account + Verif[(std::max(n_vx,0)+iloop)+2*(std::max(n_vy,0)+jloop)+4*(std::max(n_vz,0)+kloop)]=false; + } + } + } + + } + + } + } + } + + + //Loop over the future refined blocks R+1 + for (int i=0; i<2; ++i) { + for (int j=0; j<2; ++j) { + for (int k=0; k<2; ++k) { + if (Verif[i+2*j+4*k]){ + Realf Datagrossum = 0; + + for (int i2=-1; i2<2; ++i2) { + for (int j2=-1; j2<2; ++j2) { + for (int k2=-1; k2<2; ++k2) { + //Sum over the 8 cells of level R in order to reproduce the coarse cell R-1 + Datagrossum += M[1-i2*(1-2*i)][1-j2*(1-2*j)][1-k2*(1-2*k)]*Datagros[1+i+i2][1+j+j2][1+k+k2]; + } + } + } + + Realf D = abs( data[localID*WID3+cellIndex(1+i,1+j,1+k)] - Datagrossum ); + // The idea is to always compare the 8 central cells of level R with the reproduce R-1 cells + + vmesh::LocalID Indicesraf[3]; + Indicesraf[0] = 2*Indices[0]+i ; + Indicesraf[1] = 2*Indices[1]+j ; + Indicesraf[2] = 2*Indices[2]+k ; + + Realf Dcomp = cell->getVelocityBlockMinValue(popID); + + if(getObjectWrapper().particleSpecies[popID].CriteriaMethod==1){ //epsilon + Dcomp= getObjectWrapper().particleSpecies[popID].CriteriaValue; + }else if(getObjectWrapper().particleSpecies[popID].CriteriaMethod==2){ //epsilon*2^-R + Dcomp= getObjectWrapper().particleSpecies[popID].CriteriaValue/(1u << getObjectWrapper().particleSpecies[popID].RefinementLevel); + } + if (D > Dcomp){ + // We should create a new block for R+1 + int addWidthV = 1;//getObjectWrapper().particleSpecies[popID+1].sparseBlockAddWidthV; + for (int offset_vx=-addWidthV;offset_vx<=addWidthV;offset_vx++) { + for (int offset_vy=-addWidthV;offset_vy<=addWidthV;offset_vy++) { + for (int offset_vz=-addWidthV;offset_vz<=addWidthV;offset_vz++) { + const vmesh::GlobalID globalIDraf = vmeshraf->getGlobalID(Indicesraf[0]+offset_vx,Indicesraf[1]+offset_vy,Indicesraf[2]+offset_vz); + if (globalIDraf== vmeshraf->invalidGlobalID()) { + // std::cout<< " GlobalID bug not normal" << "Indices[0]+offset_vx" << Indices[0]+offset_vx << "Indices[1]+offset_vy " << Indices[1]+offset_vy << "Indices[1]+offset_vy" << Indices[2]+offset_vz <-1; --popID) { + if(getObjectWrapper().particleSpecies[popID].MaxRefinementLevel>0 && getObjectWrapper().particleSpecies[popID].RefinementLevelget_velocity_mesh(popID); + vmesh::VelocityMesh* vmeshraf = cell->get_velocity_mesh(popID+1); + + for (vmesh::GlobalID globalIDraf : ListBlockExist[popID+1]) { + + vmesh::LocalID Indicesraf[3]; + vmeshraf->getIndices(globalIDraf, Indicesraf[0], Indicesraf[1], Indicesraf[2]); + + vmesh::LocalID Indices[3]; + Indices[0] = Indicesraf[0]/2; + Indices[1] = Indicesraf[1]/2; + Indices[2] = Indicesraf[2]/2; + + int addWidthV = 1;//getObjectWrapper().particleSpecies[popID].sparseBlockAddWidthV; + for (int offset_vx=-addWidthV;offset_vx<=addWidthV;offset_vx++) { + for (int offset_vy=-addWidthV;offset_vy<=addWidthV;offset_vy++) { + for (int offset_vz=-addWidthV;offset_vz<=addWidthV;offset_vz++) { + const vmesh::GlobalID globalID = vmesh->getGlobalID(Indices[0]+offset_vx,Indices[1]+offset_vy,Indices[2]+offset_vz); + if (globalID== vmesh->invalidGlobalID()) { + std::cout<< " GlobalID bug not normal" << "Indices[0]+offset_vx" << Indices[0]+offset_vx << "Indices[1]+offset_vy " << Indices[1]+offset_vy << "Indices[1]+offset_vy" << Indices[2]+offset_vz <0){ + + vmesh::VelocityMesh* vmeshinit = cell->get_velocity_mesh(popID); + vmesh::LocalID Localsize= vmeshinit->size(); + + for (vmesh::LocalID localID=0; localIDget_velocity_mesh(popID); + vmesh::GlobalID globalID = vmeshinit2->getGlobalID(localID); + if (ListBlockExist[popID].find(globalID) == ListBlockExist[popID].end()) { + cell->remove_velocity_block(globalID,popID); + Localsize-=1; + localID-=1; + } + } + + + if(getObjectWrapper().particleSpecies[popID].RefinementLevel!=0){ + for (vmesh::GlobalID globalID : ListBlockExist[popID]) { + if(cell->add_velocity_block(globalID,popID)){ + //True if it's a new block + //need to adapt the creation to the level of refinement needed + vmesh::VelocityMesh* vmesh = cell->get_velocity_mesh(popID); + vmesh::VelocityMesh* vmeshgros = cell->get_velocity_mesh(popID-1); + Realf *datagros = cell->get_velocity_blocks(popID-1)->getData(); + Realf *data = cell->get_velocity_blocks(popID)->getData(); + + vmesh::LocalID localIDcreated=vmesh->getLocalID(globalID); + if (localIDcreated== vmesh->invalidLocalID()) { + //std::cout<< " localIDcreated bug not normal" <getIndices(globalID, Indices[0], Indices[1], Indices[2]); + // i, j and k indicate the separation of the coarse block into eight refined blocks + int i = Indices[0]%2; + int j = Indices[1]%2; + int k = Indices[2]%2; + + vmesh::LocalID Indicesgrosinit[3]; + Indicesgrosinit[0] = (Indices[0]-i)/2; + Indicesgrosinit[1] = (Indices[1]-j)/2; + Indicesgrosinit[2] = (Indices[2]-k)/2; + + for (int i2=0; i2<2; ++i2) { + for (int j2=0; j2<2; ++j2) { + for (int k2=0; k2<2; ++k2) { + //Loop over the 8 coarsed cells located on the WID3 refined cells (the refined block) + // We ensure that the cells are well initialised to 0 + data[localIDcreated*WID3+cellIndex(2*i2,2*j2,2*k2)]=0.0; + data[localIDcreated*WID3+cellIndex(2*i2+1,2*j2,2*k2)]=0.0; + data[localIDcreated*WID3+cellIndex(2*i2,2*j2+1,2*k2)]=0.0; + data[localIDcreated*WID3+cellIndex(2*i2,2*j2,2*k2+1)]=0.0; + data[localIDcreated*WID3+cellIndex(2*i2+1,2*j2+1,2*k2)]=0.0; + data[localIDcreated*WID3+cellIndex(2*i2+1,2*j2,2*k2+1)]=0.0; + data[localIDcreated*WID3+cellIndex(2*i2,2*j2+1,2*k2+1)]=0.0; + data[localIDcreated*WID3+cellIndex(2*i2+1,2*j2+1,2*k2+1)]=0.0; + for (int i4=-1; i4<2; ++i4){ + for (int j4=-1; j4<2; ++j4){ + for (int k4=-1; k4<2; ++k4){ + // Loop over the coarse neighbour cell and the actual cell + // i5, j5 and k5 are the indices of the coarse cell that will be used for interpolation + //It can be negative or >WID-1 if we need to go to another block + int i5=2*i+i2+i4; + int j5=2*j+j2+j4; + int k5=2*k+k2+k4; + + vmesh::LocalID Indicesgros[3]; + Indicesgros[0] = Indicesgrosinit[0]; + Indicesgros[1] = Indicesgrosinit[1]; + Indicesgros[2] = Indicesgrosinit[2]; + //We check whether we need to move to a different block. + if(i5<0){ + Indicesgros[0]-= 1; + }else if(i5> WID-1){ + Indicesgros[0]+= 1; + } + + if(j5<0){ + Indicesgros[1]-= 1; + }else if(j5> WID-1){ + Indicesgros[1]+= 1; + } + + if(k5<0){ + Indicesgros[2]-= 1; + }else if(k5> WID-1){ + Indicesgros[2]+= 1; + } + + vmesh::GlobalID globalIDgros=vmeshgros->getGlobalID(Indicesgros); + if (globalIDgros== vmeshgros->invalidGlobalID()) { + //std::cout<< " GlobalID bug not normal"<< "Indices[0] " << Indices[0] << "Indicesgros[0] "<< Indicesgros[0] << " GlobalID " << globalID << " PopID " << popID <getLocalID(globalIDgros); + if (localIDgros== vmeshgros->invalidLocalID()) { + //std::cout<< " localIDgros bug not normal"<< "Indices[0] " << Indices[0] << "Indicesgros[0] "<< Indicesgros[0] << " GlobalIDgros " << globalIDgros << " i5, j5, k5 " << i5 << " " << j5 << " " << k5 << " " <(WID-1) || j5<0 || j5>(WID-1) || k5<0 || k5>(WID-1) ) { + //std::cout<< " localidgros bug not normal"<< "Indices[0] " << Indices[0] << "Indicesgros[0] "<< Indicesgros[0] << " GlobalIDgros " << globalIDgros << " i5, j5, k5 " << i5 << " " << j5 << " " << k5 << " " <get_velocity_mesh(popID); + vmesh::LocalID localIDcreated=vmesh->getLocalID(globalID); + uint8_t *ghost = cell->get_velocity_blocks(popID)->getGhost(localIDcreated); + ghost[0]=1; + } + } + }else{ + for (vmesh::GlobalID globalID : ListBlockExist[popID]) { + cell->add_velocity_block(globalID,popID); + //ghost not important for this one + vmesh::VelocityMesh* vmesh = cell->get_velocity_mesh(popID); + vmesh::LocalID localIDcreated=vmesh->getLocalID(globalID); + uint8_t *ghost = cell->get_velocity_blocks(popID)->getGhost(localIDcreated); + ghost[0]=1; + } + } + + } + } +} +} + +//Every vamr_refinedStep we check all the velocity cells with the vAMR criterion with the 5th order +void RefinedOrder5(dccrg::Dccrg& mpiGrid, + const std::vector& cells){ + + + Realf M[5][5][5]; + Realf A[5] = {-3.0/128.0, 11.0/64.0, 1.0, -11.0/64.0, 3.0/128.0}; + + for (int i0 = 0; i0 < 5; i0++) { + for (int j0 = 0; j0 < 5; j0++) { + for (int k0 = 0; k0 < 5; k0++) { + M[i0][j0][k0] = A[i0] * A[j0] * A[k0]; + } + } + } + + #pragma omp parallel for schedule(dynamic,1) + for (size_t c=0; c ListBlockExist[getObjectWrapper().particleSpecies.size()]; + + for (int popID=(getObjectWrapper().particleSpecies.size()-2); popID>-1; --popID) { + if(getObjectWrapper().particleSpecies[popID].MaxRefinementLevel>0 && getObjectWrapper().particleSpecies[popID].RefinementLevelget_velocity_mesh(popID); + vmesh::VelocityMesh* vmeshraf = cell->get_velocity_mesh(popID+1); + Realf *data = cell->get_velocity_blocks(popID)->getData(); + + for (vmesh::LocalID localID=0; localIDsize(); ++localID) { + vmesh::GlobalID globalID = vmesh->getGlobalID(localID); + vmesh::LocalID Indices[3]; + vmesh->getIndices(globalID, Indices[0], Indices[1], Indices[2]); + + if (getObjectWrapper().particleSpecies[popID].RefinementLevel==0) { + ListBlockExist[popID].insert(globalID); + } + + bool Verif=true; + //We will need 2 ghost cells of level R-1 in every direction, meaning that we need a 6x6x6 grid of level R-1 cells for our block of level R + Realf Datagros[6][6][6]; + + for (int neighbour_vx=0; neighbour_vx<6 && Verif; ++neighbour_vx) { + int n_vx=neighbour_vx/2 -1; //gives the shift in block + int i=abs(neighbour_vx%2); + for (int neighbour_vy=0; neighbour_vy<6 && Verif; ++neighbour_vy) { + int n_vy=neighbour_vy/2 -1; + int j=abs(neighbour_vy%2); + for (int neighbour_vz=0; neighbour_vz<6 && Verif; ++neighbour_vz) { + int n_vz=neighbour_vz/2 -1; + int k=abs(neighbour_vz%2); + + vmesh::LocalID localIDneigh = vmesh->invalidLocalID(); + + if(n_vx==0 && n_vy==0 && n_vz==0){ + //The cell of level R-1 will be defined on the actual block of level R + localIDneigh = localID; + }else{ + //The cell of level R-1 will be defined on a neighbouring block of level R + const vmesh::GlobalID globalIDneigh = vmesh->getGlobalID(Indices[0]+n_vx,Indices[1]+n_vy,Indices[2]+n_vz); + if (globalIDneigh != vmesh->invalidGlobalID()) { + localIDneigh=vmesh->getLocalID(globalIDneigh); + } + } + + if (localIDneigh != vmesh->invalidLocalID()) { + //The block exists, meaning that we can compute the cell value for level R-1 + Realf Datagrossum=0; + for (int i2=0; i2<2; ++i2) { + for (int j2=0; j2<2; ++j2) { + for (int k2=0; k2<2; ++k2) { + Datagrossum += data[localIDneigh*WID3+cellIndex(2*i+i2,2*j+j2,2*k+k2)]; + } + } + } + //Computation of the grid R-1 + Datagros[neighbour_vx][neighbour_vy][neighbour_vz] = Datagrossum/8.0; + }else{ + //If one of the neighbouring blocks does not exist, it will directly impact all future R-1 blocks + Verif=false; + } + + } + } + } + + + if (Verif){ + for (int i=0; i<2; ++i) { + for (int j=0; j<2; ++j) { + for (int k=0; k<2; ++k) { + //Loop over the future refined blocks R+1 + Realf Datagrossum = 0; + + for (int i2=-2; i2<3; ++i2) { + for (int j2=-2; j2<3; ++j2) { + for (int k2=-2; k2<3; ++k2) { + //Sum over the 8 cells of level R in order to reproduce the coarse cell R-1 + Datagrossum += M[2-i2*(1-2*i)][2-j2*(1-2*j)][2-k2*(1-2*k)]*Datagros[2+i+i2][2+j+j2][2+k+k2]; + } + } + } + + Realf D = abs( data[localID*WID3+cellIndex(1+i,1+j,1+k)] - Datagrossum ); + // The idea is to always compare the 8 central cells of level R with the reproduce R-1 cells + vmesh::LocalID Indicesraf[3]; + Indicesraf[0] = 2*Indices[0]+i ; + Indicesraf[1] = 2*Indices[1]+j ; + Indicesraf[2] = 2*Indices[2]+k ; + + Realf Dcomp = cell->getVelocityBlockMinValue(popID); + + if(getObjectWrapper().particleSpecies[popID].CriteriaMethod==1){ //epsilon + Dcomp= getObjectWrapper().particleSpecies[popID].CriteriaValue; + }else if(getObjectWrapper().particleSpecies[popID].CriteriaMethod==2){ //epsilon*2^-R + Dcomp= getObjectWrapper().particleSpecies[popID].CriteriaValue/(1u << getObjectWrapper().particleSpecies[popID].RefinementLevel); + } + if (D > Dcomp){ + // We should create a new block for R+1 + int addWidthV = 1;//getObjectWrapper().particleSpecies[popID+1].sparseBlockAddWidthV; + for (int offset_vx=-addWidthV;offset_vx<=addWidthV;offset_vx++) { + for (int offset_vy=-addWidthV;offset_vy<=addWidthV;offset_vy++) { + for (int offset_vz=-addWidthV;offset_vz<=addWidthV;offset_vz++) { + const vmesh::GlobalID globalIDraf = vmeshraf->getGlobalID(Indicesraf[0]+offset_vx,Indicesraf[1]+offset_vy,Indicesraf[2]+offset_vz); + if (globalIDraf== vmeshraf->invalidGlobalID()) { + // std::cout<< " GlobalID bug not normal" << "Indices[0]+offset_vx" << Indices[0]+offset_vx << "Indices[1]+offset_vy " << Indices[1]+offset_vy << "Indices[1]+offset_vy" << Indices[2]+offset_vz <-1; --popID) { + if(getObjectWrapper().particleSpecies[popID].MaxRefinementLevel>0 && getObjectWrapper().particleSpecies[popID].RefinementLevelget_velocity_mesh(popID); + vmesh::VelocityMesh* vmeshraf = cell->get_velocity_mesh(popID+1); + + for (vmesh::GlobalID globalIDraf : ListBlockExist[popID+1]) { + + vmesh::LocalID Indicesraf[3]; + vmeshraf->getIndices(globalIDraf, Indicesraf[0], Indicesraf[1], Indicesraf[2]); + + vmesh::LocalID Indices[3]; + Indices[0] = Indicesraf[0]/2; + Indices[1] = Indicesraf[1]/2; + Indices[2] = Indicesraf[2]/2; + + int addWidthV = 1;//getObjectWrapper().particleSpecies[popID].sparseBlockAddWidthV; //Replace species.sparseBlockAddWidthV + for (int offset_vx=-addWidthV;offset_vx<=addWidthV;offset_vx++) { + for (int offset_vy=-addWidthV;offset_vy<=addWidthV;offset_vy++) { + for (int offset_vz=-addWidthV;offset_vz<=addWidthV;offset_vz++) { + const vmesh::GlobalID globalID = vmesh->getGlobalID(Indices[0]+offset_vx,Indices[1]+offset_vy,Indices[2]+offset_vz); + if (globalID== vmesh->invalidGlobalID()) { + std::cout<< " GlobalID bug not normal" << "Indices[0]+offset_vx" << Indices[0]+offset_vx << "Indices[1]+offset_vy " << Indices[1]+offset_vy << "Indices[1]+offset_vy" << Indices[2]+offset_vz <0){ + + vmesh::VelocityMesh* vmeshinit = cell->get_velocity_mesh(popID); + vmesh::LocalID Localsize= vmeshinit->size(); + + for (vmesh::LocalID localID=0; localIDget_velocity_mesh(popID); + vmesh::GlobalID globalID = vmeshinit2->getGlobalID(localID); + if (ListBlockExist[popID].find(globalID) == ListBlockExist[popID].end()) { + cell->remove_velocity_block(globalID,popID); + Localsize-=1; + localID-=1; + } + } + + if(getObjectWrapper().particleSpecies[popID].RefinementLevel!=0){ + for (vmesh::GlobalID globalID : ListBlockExist[popID]) { + if(cell->add_velocity_block(globalID,popID)){ + //True if it's a new block + //need to adapt the creation to the level of refinement needed + vmesh::VelocityMesh* vmesh = cell->get_velocity_mesh(popID); + vmesh::VelocityMesh* vmeshgros = cell->get_velocity_mesh(popID-1); + Realf *datagros = cell->get_velocity_blocks(popID-1)->getData(); + Realf *data = cell->get_velocity_blocks(popID)->getData(); + + vmesh::LocalID localIDcreated=vmesh->getLocalID(globalID); + if (localIDcreated== vmesh->invalidLocalID()) { + std::cout<< " localIDcreated bug not normal" <getIndices(globalID, Indices[0], Indices[1], Indices[2]); + // i, j and k indicate the separation of the coarse block into eight refined blocks + int i = Indices[0]%2; + int j = Indices[1]%2; + int k = Indices[2]%2; + + vmesh::LocalID Indicesgrosinit[3]; + Indicesgrosinit[0] = (Indices[0]-i)/2; + Indicesgrosinit[1] = (Indices[1]-j)/2; + Indicesgrosinit[2] = (Indices[2]-k)/2; + + for (int i2=0; i2<2; ++i2) { + for (int j2=0; j2<2; ++j2) { + for (int k2=0; k2<2; ++k2) { + //Loop over the 8 coarsed cells located on the WID3 refined cells (the refined block) + // We ensure that the cells are well initialised to 0 + data[localIDcreated*WID3+cellIndex(2*i2,2*j2,2*k2)]=0.0; + data[localIDcreated*WID3+cellIndex(2*i2+1,2*j2,2*k2)]=0.0; + data[localIDcreated*WID3+cellIndex(2*i2,2*j2+1,2*k2)]=0.0; + data[localIDcreated*WID3+cellIndex(2*i2,2*j2,2*k2+1)]=0.0; + data[localIDcreated*WID3+cellIndex(2*i2+1,2*j2+1,2*k2)]=0.0; + data[localIDcreated*WID3+cellIndex(2*i2+1,2*j2,2*k2+1)]=0.0; + data[localIDcreated*WID3+cellIndex(2*i2,2*j2+1,2*k2+1)]=0.0; + data[localIDcreated*WID3+cellIndex(2*i2+1,2*j2+1,2*k2+1)]=0.0; + for (int i4=-2; i4<3; ++i4){ + for (int j4=-2; j4<3; ++j4){ + for (int k4=-2; k4<3; ++k4){ + // Loop over the coarse neighbour cell and the actual cell + // i5, j5 and k5 are the indices of the coarse cell that will be used for interpolation + //It can be negative or >WID-1 if we need to go to another block + int i5=2*i+i2+i4; + int j5=2*j+j2+j4; + int k5=2*k+k2+k4; + + vmesh::LocalID Indicesgros[3]; + Indicesgros[0] = Indicesgrosinit[0]; + Indicesgros[1] = Indicesgrosinit[1]; + Indicesgros[2] = Indicesgrosinit[2]; + //We check whether we need to move to a different block + if(i5<0){ + Indicesgros[0]-= 1; + }else if(i5> WID-1){ + Indicesgros[0]+= 1; + } + + if(j5<0){ + Indicesgros[1]-= 1; + }else if(j5> WID-1){ + Indicesgros[1]+= 1; + } + + if(k5<0){ + Indicesgros[2]-= 1; + }else if(k5> WID-1){ + Indicesgros[2]+= 1; + } + + vmesh::GlobalID globalIDgros=vmeshgros->getGlobalID(Indicesgros); + if (globalIDgros== vmeshgros->invalidGlobalID()) { + std::cout<< " GlobalID bug not normal"<< "Indices[0] " << Indices[0] << "Indicesgros[0] "<< Indicesgros[0] << " GlobalID " << globalID << " PopID " << popID <getLocalID(globalIDgros); + if (localIDgros== vmeshgros->invalidLocalID()) { + std::cout<< " localIDgros bug not normal"<< "Indices[0] " << Indices[0] << "Indicesgros[0] "<< Indicesgros[0] << " GlobalIDgros " << globalIDgros << " i5, j5, k5 " << i5 << " " << j5 << " " << k5 << " " <get_velocity_mesh(popID); + vmesh::LocalID localIDcreated=vmesh->getLocalID(globalID); + uint8_t *ghost = cell->get_velocity_blocks(popID)->getGhost(localIDcreated); + ghost[0]=1; + } + } + }else{ + for (vmesh::GlobalID globalID : ListBlockExist[popID]) { + cell->add_velocity_block(globalID,popID); + //ghost not important for this one + vmesh::VelocityMesh* vmesh = cell->get_velocity_mesh(popID); + vmesh::LocalID localIDcreated=vmesh->getLocalID(globalID); + uint8_t *ghost = cell->get_velocity_blocks(popID)->getGhost(localIDcreated); + ghost[0]=1; + } + } + + } + } +} +} + +//Update the ghost parameter and remove the newly created cells that don't respect the criterion +void SmallRefinedOrder1(dccrg::Dccrg& mpiGrid, + const std::vector& cells){ + +#pragma omp parallel for schedule(dynamic,1) +for (size_t c=0; c0 && getObjectWrapper().particleSpecies[popID].RefinementLevel!=0){ + vmesh::VelocityMesh* vmesh = cell->get_velocity_mesh(popID); + Realf *data = cell->get_velocity_blocks(popID)->getData(); + uint8_t *ghost = cell->get_velocity_blocks(popID)->getGhost(); + + vmesh::LocalID Localsize= vmesh->size(); + for (vmesh::LocalID localID=0; localIDgetVelocityBlockMinValue(popID-1); + + if(getObjectWrapper().particleSpecies[popID-1].CriteriaMethod==1){ //epsilon + Dcomp= getObjectWrapper().particleSpecies[popID-1].CriteriaValue; + }else if(getObjectWrapper().particleSpecies[popID-1].CriteriaMethod==2){ //epsilon*2^-R + Dcomp= getObjectWrapper().particleSpecies[popID-1].CriteriaValue/(1u << getObjectWrapper().particleSpecies[popID-1].RefinementLevel); + } + if (D > Dcomp){ + // We keep the block + ghost[localID]=1; // This could be changed to 2 to distinguish the newly created cells from those that already exist + }else{ + //If the block don't need to exist anymore, it is removed + vmesh::VelocityMesh* vmesh = cell->get_velocity_mesh(popID); + vmesh::GlobalID globalID = vmesh->getGlobalID(localID); + cell->remove_velocity_block(globalID,popID); + Localsize-=1; + localID-=1; + } + } + } + + } + } +} +} + +//Fixed ghost to 1 for all the initial velocity cells to avoid destruction +void GhostFixation(dccrg::Dccrg& mpiGrid, + const std::vector& cells){ + +#pragma omp parallel for schedule(dynamic,1) +for (size_t c=0; cget_velocity_mesh(popID); + uint8_t *ghost = cell->get_velocity_blocks(popID)->getGhost(); + for (vmesh::LocalID localID=0; localIDsize(); ++localID) { + ghost[localID]=1; + } + } +} +} diff --git a/vlasovsolver/arch_moments.h b/vlasovsolver/arch_moments.h index c5b7c39d20..2fb06f2fbe 100644 --- a/vlasovsolver/arch_moments.h +++ b/vlasovsolver/arch_moments.h @@ -44,6 +44,11 @@ void blockVelocityFirstMoments(vmesh::VelocityBlockContainer *blockContainer, REAL (&array)[SIZE], uint nBlocks); +template +void blockVelocityFirstMomentsVamr(vmesh::VelocityBlockContainer *blockContainer, + REAL (&array)[SIZE], + uint nBlocks); + template void blockVelocitySecondMoments(vmesh::VelocityBlockContainer *blockContainer, const REAL averageVX, @@ -52,6 +57,14 @@ void blockVelocitySecondMoments(vmesh::VelocityBlockContainer *blockContainer, REAL (&array)[SIZE], uint nBlocks); +template +void blockVelocitySecondMomentsVamr(vmesh::VelocityBlockContainer *blockContainer, + const REAL averageVX, + const REAL averageVY, + const REAL averageVZ, + REAL (&array)[SIZE], + uint nBlocks); + void calculateMoments_R( dccrg::Dccrg& mpiGrid, const std::vector& cells, @@ -66,6 +79,28 @@ void calculateMoments_V( const bool initialCompute=false ); +void vamr_transfer_values( + dccrg::Dccrg& mpiGrid, + const std::vector& cells, + const uint popID); + +void changeRefined(spatial_cell::SpatialCell* cell, + const uint popID); + +void RefinedOrder1(dccrg::Dccrg& mpiGrid, + const std::vector& cells); + +void RefinedOrder3(dccrg::Dccrg& mpiGrid, + const std::vector& cells); + +void RefinedOrder5(dccrg::Dccrg& mpiGrid, + const std::vector& cells); + +void SmallRefinedOrder1(dccrg::Dccrg& mpiGrid, + const std::vector& cells); + +void GhostFixation(dccrg::Dccrg& mpiGrid, + const std::vector& cells); // ***** TEMPLATE FUNCTION DEFINITIONS ***** // @@ -107,6 +142,39 @@ void blockVelocityFirstMoments( }, array); } +/** Calculate the first velocity moments for the vAMR version **/ +template inline +void blockVelocityFirstMomentsVamr( + vmesh::VelocityBlockContainer *blockContainer, + REAL (&array)[SIZE], + uint nBlocks) { + + arch::parallel_reduce({WID, WID, WID, nBlocks}, + ARCH_LOOP_LAMBDA (const uint i, const uint j, const uint k, const uint blockLID, Real *lsum ) { + + Realf *data = blockContainer->getData(); + uint8_t *refined = blockContainer->getRefined(); + Real *blockParameters = blockContainer->getParameters(); + const Realf* avgs = &data[blockLID*WID3]; + const uint8_t* check = &refined[blockLID*WID3]; + const Real* blockParamsZ = &blockParameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS]; + const Real DV3 = blockParamsZ[BlockParams::DVX]*blockParamsZ[BlockParams::DVY]*blockParamsZ[BlockParams::DVZ]; + const Real HALF = 0.5; + + ARCH_INNER_BODY(i, j, k, blockLID, lsum) { + if (!check[cellIndex(i,j,k)]){ + const Real VX = blockParamsZ[BlockParams::VXCRD] + (i+HALF)*blockParamsZ[BlockParams::DVX]; + const Real VY = blockParamsZ[BlockParams::VYCRD] + (j+HALF)*blockParamsZ[BlockParams::DVY]; + const Real VZ = blockParamsZ[BlockParams::VZCRD] + (k+HALF)*blockParamsZ[BlockParams::DVZ]; + lsum[0] += avgs[cellIndex(i,j,k)] * DV3; + lsum[1] += avgs[cellIndex(i,j,k)]*VX * DV3; + lsum[2] += avgs[cellIndex(i,j,k)]*VY * DV3; + lsum[3] += avgs[cellIndex(i,j,k)]*VZ * DV3; + }; + }; + }, array); +} + /** Calculate the second velocity moments for the velocity blocks, and add * results to 'array', which must have at least size three. After this function * returns, the contents of 'array' are as follows: array[0]=n(Vx-Vx0); @@ -152,4 +220,42 @@ void blockVelocitySecondMoments( }, array); } +/** Calculate the second velocity moments for the vAMR version **/ +template inline +void blockVelocitySecondMomentsVamr( + vmesh::VelocityBlockContainer *blockContainer, + const REAL averageVX, + const REAL averageVY, + const REAL averageVZ, + REAL (&array)[SIZE], + uint nBlocks) { + + arch::parallel_reduce({WID, WID, WID, nBlocks}, + ARCH_LOOP_LAMBDA (const uint i, const uint j, const uint k, const uint blockLID, Real *lsum ) { + + Realf *data = blockContainer->getData(); + uint8_t *refined = blockContainer->getRefined(); + Real *blockParameters = blockContainer->getParameters(); + const Realf* avgs = &data[blockLID*WID3]; + const uint8_t* check = &refined[blockLID*WID3]; + const Real* blockParams = &blockParameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS]; + const Real DV3 = blockParams[BlockParams::DVX]*blockParams[BlockParams::DVY]*blockParams[BlockParams::DVZ]; + const Real HALF = 0.5; + + ARCH_INNER_BODY(i, j, k, blockLID, lsum) { + if (!check[cellIndex(i,j,k)]){ + const Real VX = blockParams[BlockParams::VXCRD] + (i+HALF)*blockParams[BlockParams::DVX]; + const Real VY = blockParams[BlockParams::VYCRD] + (j+HALF)*blockParams[BlockParams::DVY]; + const Real VZ = blockParams[BlockParams::VZCRD] + (k+HALF)*blockParams[BlockParams::DVZ]; + lsum[0] += avgs[cellIndex(i,j,k)] * (VX - averageVX) * (VX - averageVX) * DV3; + lsum[1] += avgs[cellIndex(i,j,k)] * (VY - averageVY) * (VY - averageVY) * DV3; + lsum[2] += avgs[cellIndex(i,j,k)] * (VZ - averageVZ) * (VZ - averageVZ) * DV3; + lsum[3] += avgs[cellIndex(i,j,k)] * (VY - averageVY) * (VZ - averageVZ) * DV3; + lsum[4] += avgs[cellIndex(i,j,k)] * (VX - averageVX) * (VZ - averageVZ) * DV3; + lsum[5] += avgs[cellIndex(i,j,k)] * (VX - averageVX) * (VY - averageVY) * DV3; + }; + }; + }, array); +} + #endif diff --git a/vlasovsolver/vlasovmover.cpp b/vlasovsolver/vlasovmover.cpp index 8ac4d5f607..e6bf414d90 100644 --- a/vlasovsolver/vlasovmover.cpp +++ b/vlasovsolver/vlasovmover.cpp @@ -342,6 +342,17 @@ void calculateSpatialTranslation( } } + if (P::activateVamr){ + phiprof::Timer vAMR_transfer {"vAMR transfer between grids in spatial"}; + for (uint popID=(getObjectWrapper().particleSpecies.size()-1); popID>0; --popID) { + if(getObjectWrapper().particleSpecies[popID].RefinementLevel>0){ + // Update communication from level R+1 to level R + vamr_transfer_values(mpiGrid,local_propagated_cells,popID-1); + } + } + vAMR_transfer.stop(); + } + if (Parameters::prepareForRebalance == true) { // clear weight on all local cells for (size_t c=0; c 0 && step < (globalMaxSubcycles - 1)) { + adjustVelocityBlocks(mpiGrid, acceleratedCells, false, popID, true); } } @@ -479,7 +485,8 @@ void calculateAcceleration(const uint popID,const uint globalMaxSubcycles,const * @param mpiGrid Parallel grid library. * @param dt Time step.*/ void calculateAcceleration(dccrg::Dccrg& mpiGrid, - const Real dt + const Real dt, + const bool ShouldRefined ) { typedef Parameters P; const vector& cells = getLocalCells(); @@ -487,18 +494,28 @@ void calculateAcceleration(dccrg::Dccrg& int myRank; MPI_Comm_rank(MPI_COMM_WORLD,&myRank); + if (P::activateVamr){ + //Update the ghost parameter and remove the newly created cells that don't respect the criterion + SmallRefinedOrder1(mpiGrid, cells); + } + if (dt == 0.0 && P::tstep > 0) { // Even if acceleration is turned off we need to adjust velocity blocks // because the boundary conditions may have altered the velocity space, // and to update changes in no-content blocks during translation. for (uint popID=0; popID& } // Accelerate population over one subcycle step calculateAcceleration(popID,(uint)globalMaxSubcycles,step,mpiGrid,acceleratedCells,dt); + if(step==0 && (uint)globalMaxSubcycles > 1){ + adjustVelocityBlocks(mpiGrid, cells, false, popID, true); + } } // for-loop over acceleration substeps // final adjust for all cells, also updating full remote block lists - adjustVelocityBlocks(mpiGrid, cells, true, popID); + adjustVelocityBlocks(mpiGrid, cells, true, popID, true); } // for-loop over particle species } + if (P::activateVamr){ + phiprof::Timer vAMR_accel_modif {"vAMR in acceleration"}; + phiprof::Timer vAMR_transfer_accel {"vAMR transfer between grids in accel"}; + for (uint popID=(getObjectWrapper().particleSpecies.size()-1); popID>0; --popID) { + if(getObjectWrapper().particleSpecies[popID].RefinementLevel>0){ + // Update communication from level R+1 to level R + vamr_transfer_values(mpiGrid,cells,popID-1); + } + } + vAMR_transfer_accel.stop(); + if (ShouldRefined){ + phiprof::Timer vAMR_refined {"vAMR refined grid"}; + //Every vamr_refinedStep we check all the velocity cells with the vAMR criterion depending on the order + RefinedOrderX(mpiGrid,cells); + for (uint popID=0; popIDparameters[CellParams::P_23_DT2] = SC->parameters[CellParams::P_23]; } // for-loop over spatial cells } + +void RefinedOrderX(dccrg::Dccrg& mpiGrid, + const std::vector& cells){ + if (P::vAMRorder==1){ + RefinedOrder1(mpiGrid,cells); + }else if(P::vAMRorder==3){ + RefinedOrder3(mpiGrid,cells); + }else if(P::vAMRorder==5){ + RefinedOrder5(mpiGrid,cells); + }else { + RefinedOrder1(mpiGrid,cells); + } +} + +void FixedGhost(dccrg::Dccrg& mpiGrid, + const std::vector& cells){ + GhostFixation(mpiGrid,cells); +} diff --git a/vlasovsolver/vlasovmover.h b/vlasovsolver/vlasovmover.h index 281ededa17..6b137f4967 100644 --- a/vlasovsolver/vlasovmover.h +++ b/vlasovsolver/vlasovmover.h @@ -35,7 +35,8 @@ using namespace spatial_cell; void calculateAcceleration( dccrg::Dccrg& mpiGrid, - const Real dt + const Real dt, + const bool ShouldRefined=false ); void calculateSpatialTranslation( @@ -83,7 +84,11 @@ void calculateInterpolatedVelocityMoments( */ void calculateInitialVelocityMoments(dccrg::Dccrg& mpiGrid); +void RefinedOrderX(dccrg::Dccrg& mpiGrid, + const std::vector& cells); +void FixedGhost(dccrg::Dccrg& mpiGrid, + const std::vector& cells); #endif