Merge pull request #519 from alhom/eVlasiator

Proportionally increased substepping for electron solver around f_pe = f_ce, with a hard upper limit; not used to lower substepping so far.

parameters.cpp

@@ -111,6 +111,7 @@ Real P::maxWaveVelocity = 0.0;
 uint P::maxFieldSolverSubcycles = 0.0;
 int P::maxSlAccelerationSubcycles = 0.0;
 bool P::ResolvePlasmaPeriod = false;
+int P::maxResonantSubcycleFactor = 100;
 Real P::resistivity = NAN;
 bool P::fieldSolverDiffusiveEterms = true;
 uint P::ohmHallTerm = 0;
@@ -228,6 +229,8 @@ bool Parameters::addParameters(){
    Readparameters::add("vlasovsolver.maxCFL","The maximum CFL limit for vlasov propagation in ordinary space. Used to set timestep if dynamic_timestep is true.",0.99);
    Readparameters::add("vlasovsolver.minCFL","The minimum CFL limit for vlasov propagation in ordinary space. Used to set timestep if dynamic_timestep is true.",0.8);
    Readparameters::add("vlasovsolver.ResolvePlasmaPeriod","Require semi-lagrangian solver to resolve plasma oscillation (default false)",false);
+   Readparameters::add("vlasovsolver.maxResonantSubcycleFactor","Resolve (electron) plasma oscillation-cyclotron resonance with additional subcycles, max multiplicative substep number (default 100)",100);
+
 
    // Load balancing parameters
    Readparameters::add("loadBalance.algorithm", "Load balancing algorithm to be used", string("RCB"));
@@ -607,6 +610,7 @@ bool Parameters::getParameters(){
    Readparameters::get("vlasovsolver.maxCFL",P::vlasovSolverMaxCFL);
    Readparameters::get("vlasovsolver.minCFL",P::vlasovSolverMinCFL);
    Readparameters::get("vlasovsolver.ResolvePlasmaPeriod",P::ResolvePlasmaPeriod);
+   Readparameters::get("vlasovsolver.maxResonantSubcycleFactor",P::maxResonantSubcycleFactor);
 
 
    // Get load balance parameters

parameters.h

@@ -112,6 +112,7 @@ struct Parameters {
    static Real maxSlAccelerationRotation; /*!< Maximum rotation in acceleration for semilagrangian solver*/
    static int maxSlAccelerationSubcycles; /*!< Maximum number of subcycles in acceleration*/
    static bool ResolvePlasmaPeriod; /*!< Should the plasma period be considered in calculating dt*/
+   static int maxResonantSubcycleFactor; /*!< For cells where plasma freq~cyclotron freq, increase substepping with maximum of this factor*/
 
    static Real hallMinimumRhom;  /*!< Minimum mass density value used in the field solver.*/
    static Real hallMinimumRhoq;  /*!< Minimum charge density value used for the Hall and electron pressure gradient terms in the Lorentz force and in the field solver.*/

vlasovsolver/cpu_acc_transform.cpp

@@ -141,6 +141,17 @@ Eigen::Transform<Real,3,Eigen::Affine> compute_acceleration_transformation(
       unsigned int transformation_substeps_2; 
       transformation_substeps_2 = fabs(dt) / fabs(plasma_period*(0.1/360.0));
       transformation_substeps = transformation_substeps_2 > transformation_substeps ? transformation_substeps_2 : transformation_substeps;
+
+      // Account for extra substeps when plasma period ~ gyro period
+      Real logFactor = 1.0/fabs(log( fabs(gyro_period)/fabs(plasma_period) ));
+      if(std::isnan(logFactor) == false){
+         //constraint: max n times more subcycles
+         logFactor = min(logFactor, (Real)P::maxResonantSubcycleFactor);
+         //constraint: do not lower subcycle count
+         logFactor = max(logFactor, 1.0);
+         cout << logFactor << "\n";
+         transformation_substeps = int(transformation_substeps * logFactor);
+      }
    }
    if ((transformation_substeps < 1) && (fabs(dt)>0)) transformation_substeps=1;
 
@@ -206,9 +217,9 @@ Eigen::Transform<Real,3,Eigen::Affine> compute_acceleration_transformation(
 	 isn't needed.
       */
       if (!smallparticle) {
-	 rotation_pivot[0]-= hallPrefactor*(dBZdy - dBYdz);
-	 rotation_pivot[1]-= hallPrefactor*(dBXdz - dBZdx);
-	 rotation_pivot[2]-= hallPrefactor*(dBYdx - dBXdy);
+         rotation_pivot[0]-= hallPrefactor*(dBZdy - dBYdz);
+         rotation_pivot[1]-= hallPrefactor*(dBXdz - dBZdx);
+         rotation_pivot[2]-= hallPrefactor*(dBYdx - dBXdy);
       }
 
       // Calculate EJE only for the electron population
@@ -238,7 +249,7 @@ Eigen::Transform<Real,3,Eigen::Affine> compute_acceleration_transformation(
       } // end if (smallparticle==true) and dt>0
 
       if (smallparticle) {	  
-	  total_transform=Translation<Real,3>(deltaV) * total_transform;
+         total_transform=Translation<Real,3>(deltaV) * total_transform;
       }
 
       /* Evaluate electron pressure gradient term. This is treated as a simple