fmihpc · MaximeII · Mar 4, 2021 · Mar 5, 2021 · Mar 10, 2021 · Apr 15, 2021
diff --git a/Makefile b/Makefile
@@ -155,12 +155,12 @@ OBJS = 	version.o memoryallocation.o backgroundfield.o quadr.o dipole.o linedipo
 	sysboundary.o sysboundarycondition.o particle_species.o\
 	project.o projectTriAxisSearch.o read_gaussian_population.o\
 	Alfven.o Diffusion.o Dispersion.o Distributions.o Firehose.o\
-	Flowthrough.o Fluctuations.o Harris.o KHB.o Larmor.o Magnetosphere.o MultiPeak.o\
-	VelocityBox.o Riemann1.o Shock.o Template.o test_fp.o testHall.o test_trans.o\
+	Flowthrough.o Fluctuations.o Harris.o KHB.o Larmor.o Magnetosphere.o MultiPeak.o LossCone.o\
+	VelocityBox.o velocity_space_diffusion.o Riemann1.o Shock.o Template.o test_fp.o testHall.o test_trans.o\
 	IPShock.o object_wrapper.o\
 	verificationLarmor.o Shocktest.o grid.o ioread.o iowrite.o vlasiator.o logger.o\
 	common.o parameters.o readparameters.o spatial_cell.o\
-	vlasovmover.o $(FIELDSOLVER).o fs_common.o fs_limiters.o gridGlue.o
+	vlasovmover.o $(FIELDSOLVER).o fs_common.o fs_limiters.o gridGlue.o 
 
 # Include autogenerated dependency files, if they exist
 -include $(OBJS:%.o=%.d)
@@ -226,7 +226,7 @@ version.cpp: FORCE
 # for all files in the datareduction/ dir
 %.o: datareduction/%.cpp ${DEPS_COMMON} datareduction/datareductionoperator.h fieldtracing/fieldtracing.h sysboundary/ionosphere.h datareduction/dro_populations.h
 	@echo [CC] $<
-	$(SILENT)${CMP} ${CXXFLAGS} ${FLAGS} ${MATHFLAGS} -c $< ${INC_DCCRG} ${INC_ZOLTAN} ${INC_MPI} ${INC_BOOST} ${INC_EIGEN} ${INC_VLSV} ${INC_FSGRID}
+	$(SILENT)${CMP} ${CXXFLAGS} ${FLAGS} ${MATHFLAGS} -c $< ${INC_DCCRG} ${INC_ZOLTAN} ${INC_MPI} ${INC_BOOST} ${INC_EIGEN} ${INC_VLSV} ${INC_FSGRID} ${INC_VECTORCLASS}
 
 # for all files in the sysboundary/ dir
 %.o: sysboundary/%.cpp ${DEPS_COMMON} sysboundary/%.h backgroundfield/backgroundfield.h projects/project.h fieldsolver/fs_limiters.h

diff --git a/NU0BOX.DAT b/NU0BOX.DAT
diff --git a/common.h b/common.h
@@ -174,6 +174,18 @@ namespace CellParams {
       P_11_V,   /*!< P_xx component after propagation in velocity space */
       P_22_V,   /*!< P_yy component after propagation in velocity space */
       P_33_V,   /*!< P_zz component after propagation in velocity space */
+      P_12,     /*!< Pressure P_xy component, computed by Vlasov propagator. */
+      P_13,     /*!< Pressure P_xz component, computed by Vlasov propagator. */
+      P_23,     /*!< Pressure P_yz component, computed by Vlasov propagator. */
+      P_12_DT2, /*!< Intermediate step value for RK2 time stepping in field solver. Computed from P_12_R and P_12_V. */
+      P_13_DT2, /*!< Intermediate step value for RK2 time stepping in field solver. Computed from P_13_R and P_13_V. */
+      P_23_DT2, /*!< Intermediate step value for RK2 time stepping in field solver. Computed from P_23_R and P_23_V. */
+      P_12_R,   /*!< P_xy component after propagation in ordinary space */
+      P_13_R,   /*!< P_xz component after propagation in ordinary space */
+      P_23_R,   /*!< P_yz component after propagation in ordinary space */
+      P_12_V,   /*!< P_xy component after propagation in velocity space */
+      P_13_V,   /*!< P_xz component after propagation in velocity space */
+      P_23_V,   /*!< P_yz component after propagation in velocity space */
       EXVOL,    /*!< Volume electric field averaged over spatial cell, x-component.*/
       EYVOL,    /*!< Volume electric field averaged over spatial cell, y-component.*/
       EZVOL,    /*!< Volume electric field averaged over spatial cell, z-component.*/
@@ -214,6 +226,7 @@ namespace CellParams {
       BULKV_FORCING_Y, /*! Externally forced drift velocity (ex. from the ionosphere) */
       BULKV_FORCING_Z, /*! Externally forced drift velocity (ex. from the ionosphere) */
       FORCING_CELL_NUM, /*! Number of boundary cells that have forced a bulkv here */
+      NU0, /*!< nu0 value for subgrid diffusion */
       N_SPATIAL_CELL_PARAMS
    };
 }
@@ -475,7 +488,9 @@ RK_ORDER2_STEP1,   /*!< Two-step second order method, first step */
 RK_ORDER2_STEP2    /*!< Two-step second order method, second step */
 };
 
-const int WID = 4;         /*!< Number of cells per coordinate in a velocity block. Only a value of 4 supported by vectorized Leveque solver */
+#ifndef WID
+#define WID (4)            /*!< Number of cells per coordinate in a velocity block. Defaults to the historical 4. */
+#endif
 const int WID2 = WID*WID;  /*!< Number of cells per 2D slab in a velocity block. */
 const int WID3 = WID2*WID; /*!< Number of cells in a velocity block. */
 

diff --git a/datareduction/datareducer.cpp b/datareduction/datareducer.cpp
@@ -368,6 +368,11 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti
             continue;
          }
       }
+      if(P::systemWriteAllDROs || lowercase == "vg_nu0" || lowercase == "nu0") { // nu0 for sub-grid diffusion
+         outputReducer->addOperator(new DRO::DataReductionOperatorCellParams("vg_nu0",CellParams::NU0,1));
+         outputReducer->addMetadata(outputReducer->size()-1,"1/s","$\\mathrm{s}^{-1}$","$\\nu_0$","1.0");
+         continue;
+      }
       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];
@@ -513,6 +518,16 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti
             continue;
          }
       }
+      if(P::systemWriteAllDROs || lowercase == "populations_1dmuspace" || lowercase == "populations_vg_1dmuspace") {
+         // 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;
+            outputReducer->addOperator(new DRO::VariableMuSpace(i));
+            outputReducer->addMetadata(outputReducer->size()-1,"1/m^3","$\\mathrm{m}^{-3}$","$f(\\mu)$","1.0");
+         }
+         continue;
+      }
       if(P::systemWriteAllDROs || lowercase == "populations_precipitationlineflux" || lowercase == "populations_vg_precipitationlinedifferentialflux" || lowercase == "populations_precipitationlinedifferentialflux") {
          // Per-population precipitation differential flux (along line)
          for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) {

diff --git a/datareduction/datareductionoperator.cpp b/datareduction/datareductionoperator.cpp
@@ -20,6 +20,7 @@
  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
  */
 
+#include "vectorclass.h"
 #include <cstdlib>
 #include <mpi.h>
 #include <iostream>
@@ -28,6 +29,7 @@
 #include "datareductionoperator.h"
 #include "../object_wrapper.h"
 
+
 using namespace std;
 
 typedef Parameters P;
@@ -1651,6 +1653,97 @@ namespace DRO {
       return true;
    }
 
+   /*! \brief V-space flatten into 1D mu distribution
+    */
+   VariableMuSpace::VariableMuSpace(cuint _popID): DataReductionOperatorHasParameters(),popID(_popID) {popName = getObjectWrapper().particleSpecies[popID].name;}
+   VariableMuSpace::~VariableMuSpace() { }
+
+   std::string VariableMuSpace::getName() const {return popName + "/vg_1dmuspace";}
+
+   bool VariableMuSpace::getDataVectorInfo(std::string& dataType,unsigned int& dataSize,unsigned int& vectorSize) const {
+      dataType = "float";
+      dataSize =  sizeof(Real);
+      vectorSize = Parameters::PADmubins; //Number of bins to build muSpace
+      return true;
+   }
+
+   bool VariableMuSpace::reduceData(const SpatialCell* cell,char* buffer) {
+
+
+       const Real* parameters = cell->get_block_parameters(popID);
+       std::vector<Real>  fmu   (Parameters::PADmubins);
+       std::vector<int>   fcount(Parameters::PADmubins);
+
+       std::array<Real,3> bulkV = {cell->parameters[CellParams::VX], cell->parameters[CellParams::VY], cell->parameters[CellParams::VZ]};
+
+       std::array<Real,3> B = {cell->parameters[CellParams::PERBXVOL] +  cell->parameters[CellParams::BGBXVOL],
+                                 cell->parameters[CellParams::PERBYVOL] +  cell->parameters[CellParams::BGBYVOL],
+	                         cell->parameters[CellParams::PERBZVOL] +  cell->parameters[CellParams::BGBZVOL]};
+
+       Real Bnorm           = sqrt(B[0]*B[0] + B[1]*B[1] + B[2]*B[2]);
+       std::array<Real,3> b = {B[0]/Bnorm, B[1]/Bnorm, B[2]/Bnorm};
+
+       // Build 2d array of f(v,mu)
+       for (vmesh::LocalID n=0; n<cell->get_number_of_velocity_blocks(popID); n++) { // Iterate through velocity blocks
+	   for (uint k = 0; k < WID; ++k) for (uint j = 0; j < WID; ++j) for (uint i = 0; i < WID; ++i) { // Iterate through coordinates (z,y,x)
+
+	       //Get velocity space coordinates                    
+	       const Real VX(parameters[n * BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::VXCRD] 
+		       + (i + 0.5)*parameters[n * BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::DVX]);
+
+	       const Real VY(parameters[n * BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::VYCRD] 
+		       + (j + 0.5)*parameters[n * BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::DVY]);
+
+	       const Real VZ(parameters[n * BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::VZCRD]
+		       + (k + 0.5)*parameters[n * BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::DVZ]);
+
+	       std::array<Real,3> V = {VX,VY,VZ}; // Velocity in the cell, in the simulation frame
+
+	       std::array<Real,3> Vplasma; // Velocity in the cell, in the plasma frame
+
+	       for (int indx = 0; indx < 3; indx++) { Vplasma[indx] = (V[indx] - bulkV[indx]); }
+
+	       Real normV = sqrt(Vplasma[0]*Vplasma[0] + Vplasma[1]*Vplasma[1] + Vplasma[2]*Vplasma[2]);
+
+	       Real Vpara = Vplasma[0]*b[0] + Vplasma[1]*b[1] + Vplasma[2]*b[2];
+
+	       Real mu = Vpara/(normV+std::numeric_limits<Realf>::min()); // + min value to avoid division by 0
+
+	       int muindex;
+               Realf dmubins = 2.0 / Parameters::PADmubins; 
+	       muindex = floor((mu+1.0) / dmubins); 
+
+	       const Realf CellValue = cell->get_data(n,popID)[i+WID*j+WID*WID*k];
+
+               const Real DVX = parameters[n * BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::DVX];
+               const Real DVY = parameters[n * BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::DVY];
+               const Real DVZ = parameters[n * BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::DVZ];
+
+	       fmu   [muindex] += CellValue * DVX*DVY*DVZ / dmubins;
+
+
+	   } // End coordinates
+       } // End blocks
+
+       // Divide f by count (independent of v but needs to be computed for all mu before derivatives)
+       //for(int indmu = 0; indmu < Parameters::PADmubins; indmu++) { 
+       //    if (fcount[indmu] == 0 || fmu[indmu] <= 0.0) { fmu[indmu] = std::numeric_limits<Realf>::min();}
+       //    else {fmu[indmu] = fmu[indmu] / fcount[indmu];} 
+       //}
+
+      const char* ptr = reinterpret_cast<const char*>(fmu.data());
+      for (uint i = 0; i < Parameters::PADmubins*sizeof(Real); ++i) buffer[i] = ptr[i];
+      return true;
+   }
+
+   bool VariableMuSpace::setSpatialCell(const SpatialCell* cell) {
+      return true;
+   }
+
+   bool VariableMuSpace::writeParameters(vlsv::Writer& vlsvWriter) {
+       return true;
+   }
+
    /*! \brief Precipitation directional differential number flux (along line)
     * Evaluation of the precipitating differential flux (per population) targeted at low energies.
     * In a selected number (default: 16) of logarithmically spaced energy bins, the average of

diff --git a/datareduction/datareductionoperator.h b/datareduction/datareductionoperator.h
@@ -614,6 +614,24 @@ namespace DRO {
       std::vector<Real> channels, dataDiffFlux;
    };
 
+
+   // V-space flatten into 1D mu-distribution
+   class VariableMuSpace: public DataReductionOperatorHasParameters {
+   public:
+      VariableMuSpace(cuint popID);
+      virtual ~VariableMuSpace();
+
+      virtual bool getDataVectorInfo(std::string& dataType,unsigned int& dataSize,unsigned int& vectorSize) const;
+      virtual std::string getName() const;
+      virtual bool reduceData(const SpatialCell* cell,char* buffer);
+      virtual bool setSpatialCell(const SpatialCell* cell);
+      virtual bool writeParameters(vlsv::Writer& vlsvWriter);
+
+   protected:
+      uint popID;
+      std::string popName;
+   };
+
    // Precipitation directional differential number flux (along line)
    class VariablePrecipitationLineDiffFlux: public DataReductionOperatorHasParameters {
    public:
@@ -625,7 +643,6 @@ namespace DRO {
       virtual bool reduceData(const SpatialCell* cell,char* buffer);
       virtual bool setSpatialCell(const SpatialCell* cell);
       virtual bool writeParameters(vlsv::Writer& vlsvWriter);
-
    protected:
       uint popID;
       std::string popName;

diff --git a/item_storage.h b/item_storage.h
@@ -25,7 +25,6 @@
 
 #include <map>
 #include <string>
-
 #include "definitions.h"
 
 /** A generic storage class for storing items, such as variable values or 

diff --git a/object_wrapper.cpp b/object_wrapper.cpp
@@ -1,6 +1,7 @@
 #include "object_wrapper.h"
 #include "readparameters.h"
 #include <iostream>
+#include <string>
 
 bool ObjectWrapper::addParameters() {
    typedef Readparameters RP;

diff --git a/object_wrapper.h b/object_wrapper.h
@@ -32,6 +32,7 @@
 #include "particle_species.h"
 #include "projects/project.h"
 #include "velocity_mesh_parameters.h"
+#include <string>
 #include "sysboundary/sysboundary.h"
 
 struct ObjectWrapper {

diff --git a/parameters.cpp b/parameters.cpp
@@ -63,7 +63,8 @@ uint P::zcells_ini = numeric_limits<uint>::max();
 
 Real P::t = 0;
 Real P::t_min = 0;
-Real P::t_max = LARGE_REAL;
+Real P::t_max  = LARGE_REAL;
+Real P::dt_ceil = -1.0; 
 Real P::dt = NAN;
 Real P::vlasovSolverMaxCFL = NAN;
 Real P::vlasovSolverMinCFL = NAN;
@@ -186,7 +187,14 @@ std::vector<Realf> P::amrBoxCenterY;
 std::vector<Realf> P::amrBoxCenterZ;
 std::vector<int> P::amrBoxMaxLevel;
 vector<string> P::blurPassString;
-std::vector<int> P::numPasses; //numpasses
+vector<int> P::numPasses;
+bool P::artificialPADiff;
+Realf P::PADcoefficient;
+Realf P::PADCFL;
+int P::PADvbins;
+int P::PADmubins;
+string P::PADnu0 = string("");
+Realf P::PADfudge;
 
 std::array<FsGridTools::Task_t,3> P::manualFsGridDecomposition = {0,0,0};
 std::array<FsGridTools::Task_t,3> P::overrideReadFsGridDecomposition = {0,0,0};
@@ -313,6 +321,9 @@ bool P::addParameters() {
    RP::add("gridbuilder.timestep_max",
            "Max. value for timesteps. If t_max limit is hit first, this step will never be reached",
            numeric_limits<uint>::max());
+   RP::add("gridbuilder.dt_ceil",
+           "Maximum simulation dt in seconds.",
+           -1.0);
 
    // Field solver parameters
    RP::add("fieldsolver.maxWaveVelocity",
@@ -497,6 +508,15 @@ bool P::addParameters() {
    RP::add("AMR.transShortPencils", "if true, use one-cell pencils", false);
    RP::addComposing("AMR.filterpasses", string("AMR filter passes for each individual refinement level"));
 
+   // Diffusion parameters
+   RP::add("PAD.enable","Enable Artificial pitch-angle diffusion",0);
+   RP::add("PAD.coefficient","Set artificial pitch-angle diffusion coefficient",0);
+   RP::add("PAD.CFL","Set CFL condition",0);
+   RP::add("PAD.vbins","number of bins for velocity",200);
+   RP::add("PAD.mubins","number of bins for mu",30);
+   RP::add("PAD.file","Path of txt file for nu0", string("NU0BOX.DAT"));
+   RP::add("PAD.fudge","fudge factor to lower diffusion",0.25);
+   // Fieldtracing
    RP::add("fieldtracing.fieldLineTracer", "Field line tracing method to use for coupling ionosphere and magnetosphere (options are: Euler, BS)", std::string("Euler"));
    RP::add("fieldtracing.tracer_max_allowed_error", "Maximum allowed error for the adaptive field line tracers ", 1000);
    RP::add("fieldtracing.tracer_max_attempts", "Maximum allowed attempts for the adaptive field line tracers", 100);
@@ -878,6 +898,8 @@ void Parameters::getParameters() {
    RP::get("gridbuilder.t_max", P::t_max);
    RP::get("gridbuilder.timestep_max", P::tstep_max);
 
+   RP::get("gridbuilder.dt_ceil", P::dt_ceil);
+
    if (P::dynamicTimestep)
       P::dt = 0.0; // if dynamic timestep then first dt is always 0
 
@@ -970,6 +992,14 @@ void Parameters::getParameters() {
       P::systemWrites.push_back(0);
    }
 
+   RP::get("PAD.enable", P::artificialPADiff);   
+   RP::get("PAD.coefficient", P::PADcoefficient);
+   RP::get("PAD.CFL",P::PADCFL);
+   RP::get("PAD.vbins",P::PADvbins);
+   RP::get("PAD.mubins",P::PADmubins);
+   RP::get("PAD.file",P::PADnu0);
+   RP::get("PAD.fudge",P::PADfudge);
+
    RP::get("fieldtracing.fieldLineTracer", tracerString);
    RP::get("fieldtracing.tracer_max_allowed_error", FieldTracing::fieldTracingParameters.max_allowed_error);
    RP::get("fieldtracing.tracer_max_attempts", FieldTracing::fieldTracingParameters.max_field_tracer_attempts);
@@ -994,3 +1024,5 @@ void Parameters::getParameters() {
       abort();
    }
 }
+
+
diff --git a/parameters.h b/parameters.h
@@ -56,6 +56,8 @@ struct Parameters {
    static Real t_max; /*!< Maximum simulation time. */
    static Real dt;    /*!< The value of the timestep to use in propagation. If CflLimit defined then it is dynamically
                          updated during simulation*/
+   static Real dt_ceil; /*!< The maximum value of the timestep to use in propagation. */
+
    static Real vlasovSolverMaxCFL;   /*!< The maximum CFL limit for propagation of distribution function. Used to set
                                         timestep if useCFLlimit is true. */
    static Real vlasovSolverMinCFL;   /*!< The minimum CFL limit for propagation of distribution function. Used to set
@@ -218,6 +220,13 @@ struct Parameters {
    static bool amrTransShortPencils;        /*!< Use short or longpencils in AMR translation.*/
    static std::vector<std::string> blurPassString;
    static std::vector<int> numPasses;
+   static bool artificialPADiff; // Artificial velocity space Diffusion
+   static Realf PADcoefficient; // Artificial pitch-angle diffusion coefficient
+   static Realf PADCFL; // Artificial pitch-angle diffusion CFL
+   static int PADvbins; // Number of bins in velocity for pitch-angle diffusion
+   static int PADmubins; // Number of bins in mu for pitch-angle diffusion
+   static std::string PADnu0; // Path to txt file for nu0
+   static Realf PADfudge; // Fudge factore for diffusion   
 
    static std::array<FsGridTools::Task_t,3> manualFsGridDecomposition;
    static std::array<FsGridTools::Task_t,3> overrideReadFsGridDecomposition;