Fix diffusion for groups of particles

Author: contuandrea

include/deposits/deposit.h

@@ -38,7 +38,7 @@ namespace deposit{
         std::string _cached_path;
         RunningStateHost_t _cached_data;
         size_t _cached_nparticles;
-        size_t cached_group;
+        size_t _cached_group;
         double _cached_length;
 
         std::string _path;
@@ -62,13 +62,13 @@ namespace deposit{
         //initialise
 
         void init(singleconf& cf){
-            this->init(std::get<std::string>(cf.get("path")), std::get<size_t>(cf.get("nparticles")),std::get<double>(cf.get("length"),,std::get<size_t>(cf.get("group")));
+            this->init(std::get<std::string>(cf.get("path")), std::get<size_t>(cf.get("nparticles")),std::get<double>(cf.get("length")),std::get<size_t>(cf.get("group")));
         }
 
         void init(std::string path, size_t nparticles, double length, size_t group=DEFAULT_GROUP){
             if(_path==DEFAULT_PATH){
                 if(nparticles!=_cached_nparticles || length!=_cached_length || group!=_cached_group){
-                    _data = this->gendummy(nparticles,length);
+                    this->gendummy(nparticles,length,group);
                 }
                 else{
                     _data.resize(_nparticles);
@@ -110,11 +110,10 @@ namespace deposit{
             hydra::Random<> Generator( std::chrono::system_clock::now().time_since_epoch().count());
             size_t np=(nparticles<=0) ? MAXPARTICLES : nparticles;
             if(np==1){
-                auto data_d1=RunningStateHost_t(np,RunningTuple_t(std::copysign(1.,length),0.,0.,0., 0.,1.,0.,0.,0.,0.,0.,0.));
-                return data_d1;
+                _data=RunningStateHost_t(np,RunningTuple_t(std::copysign(1.,length),0.,0.,0., 0.,1.,0.,0.,0.,0.,0.,0.));
             }
             size_t halfsize=np/2;
-            RunningStateHost_t data_d(halfsize*2);
+            _data.resize(halfsize*2);
             auto data_de=RunningStateHost_t(halfsize,RunningTuple_t(-1.,0.,0.,0., 1.,0.,0.,0.,0.,0.,0.,0.));
             auto data_dh=RunningStateHost_t(halfsize,RunningTuple_t(1.,0.,0.,0., 1.,0.,0.,0.,0.,0.,0.,0.));
 
@@ -127,17 +126,18 @@ namespace deposit{
                 Generator.SetSeed(seed);
                 Generator.Uniform(0., length, data_dh.begin(_tc_z), data_dh.end(_tc_z));
             }
-            hydra::copy(data_de,hydra::make_range(data_d.begin(),data_d.begin()+halfsize));
-            hydra::copy(data_dh,hydra::make_range(data_d.begin()+halfsize,data_d.end()));
+            hydra::copy(data_de,hydra::make_range(_data.begin(),_data.begin()+halfsize));
+            hydra::copy(data_dh,hydra::make_range(_data.begin()+halfsize,_data.end()));
+            
         }
 
         //get deposit
         RunningStateHost_t getdata(){return _data;}
 
         //get deposit and change if different
-        void SendToDevice(RunningStateDev_t& todev, singleconf)
-
-        }
+//         void SendToDevice(RunningStateDev_t& todev, singleconf)
+// 
+//         }
     };
 }
 

include/simulation/Evolve.h

@@ -226,7 +226,9 @@ namespace evolve{
             hydra::Vector3R wfvec(hydra::get<5>(entry8),hydra::get<6>(entry8),hydra::get<7>(entry8));
 
             //calculate diffusion velocity
-            double s = sqrt(2*(_fDiff*mobility)/_fStep)*hydra::get<_tc_gauss_x>(p);
+            //
+            //
+            double s = sqrt(2*(_fDiff*mobility)/_fStep/fabs(hydra::get<_tc_charge>(p)))*hydra::get<_tc_gauss_x>(p); //divide by square root of group size
 
             hydra::Vector3R ref(s*sin(hydra::get<_tc_angle_2>(p))*cos(hydra::get<_tc_angle_1>(p)),s*sin(hydra::get<_tc_angle_2>(p))*sin(hydra::get<_tc_angle_1>(p)),s*cos(hydra::get<_tc_angle_2>(p)));
 
@@ -358,7 +360,7 @@ namespace evolve{
             auto wfvec = interpolITvec(vvec, _beginmap_wf, ind8);
 
             //calculate diffusion velocity
-            double s = sqrt(2*(_fDiff*mobility)/_fStep)*hydra::get<_tc_gauss_x>(p);
+            double s = sqrt(2*(_fDiff*mobility)/_fStep/fabs(hydra::get<_tc_charge>(p)))*hydra::get<_tc_gauss_x>(p); //divide by square root of group size
 
             hydra::Vector3R ref(s*sin(hydra::get<_tc_angle_2>(p))*cos(hydra::get<_tc_angle_1>(p)),s*sin(hydra::get<_tc_angle_2>(p))*sin(hydra::get<_tc_angle_1>(p)),s*cos(hydra::get<_tc_angle_2>(p)));