Fix slowness in cuda backend

Author: contuandrea

include/details/datatypes.h

@@ -53,9 +53,15 @@ using ReducedTuple_t        =   hydra::tuple<double,double,double,double,double,
 ReducedTuple_t ReducedTuple_init(0.,0.,0.,0.,0.,0.,0.,0.);
 
 template<typename T>
-using VecHost_t             =   hydra::host::vector<T>; //vector container double
+using VecHost_t             =   hydra::host::vector<T>; //vector container
 template<typename T>
-using VecDev_t              =   hydra::device::vector<T>; //vector container double
+using VecDev_t              =   hydra::device::vector<T>; //vector container
+
+//maps multiarrays
+using PhysMapHost_t         =   hydra::multiarray<double,8,hydra::host::sys_t>; //map container
+using VectorMapHost_t         =   hydra::multiarray<double,3,hydra::host::sys_t>; //map container
+using PhysMapDev_t         =   hydra::multiarray<double,8,hydra::device::sys_t>; //map container
+using VectorMapDev_t         =   hydra::multiarray<double,3,hydra::device::sys_t>; //map container
 
 //currents vector
 using CurrentVec_t          =   std::vector<ReducedTuple_t>;

include/maps/preparemaps.h

@@ -434,6 +434,7 @@ namespace maps{
     class mapconfig{
     private:
         bool _multimap = false;
+        bool _isloaded = false;
         std::tuple<double,double,double,double,double,double> _borders;
 
         maps::physmap<double,11> _pm;
@@ -443,12 +444,15 @@ namespace maps{
         maps::physmap<double,4> _hmobm;
 
         VecDev_t<double> _xvec, _yvec, _zvec;
+        PhysMapHost_t _vec_map;
+        
         //efield vecs
         VecDev_t<double> _xvec_ef, _yvec_ef, _zvec_ef;
-        VecDev_t<double> efieldvecx, efieldvecy, efieldvecz;
+        VectorMapHost_t _vec_ef;
+        
         //wfield vecs
         VecDev_t<double> _xvec_wf, _yvec_wf, _zvec_wf;
-        VecDev_t<double> wfieldvecx, wfieldvecy, wfieldvecz;
+        VectorMapHost_t _vec_wf;
 
         //emob vecs
         VecDev_t<double> _xvec_emob, _yvec_emob, _zvec_emob;
@@ -459,14 +463,23 @@ namespace maps{
         VecDev_t<double> hmobvec;
 
     public:
-        mapconfig(){}
+        mapconfig() {
+            _isloaded = false;
+        };
+        
+        mapconfig(std::string map_path) {
+            this->load(map_path);
+        }
+        mapconfig(std::string efield_path, std::string emob_path, std::string hmob_path, std::string wfield_path){
+            this->load(efield_path, emob_path, hmob_path, wfield_path);
+        }
 
         void load(std::string map_path){
             INFO_LINE("Loading single map")
             _pm.init(map_path);
             _pm.prepare();
             _pm.fillspacepoints<VecDev_t<double> >(_xvec,_yvec,_zvec);
-            _pm.fillallquantities<VecDev_t<double> >(efieldvecx,efieldvecy,efieldvecz,emobvec,hmobvec,wfieldvecx,wfieldvecy,wfieldvecz);
+            _pm.fillallquantities<PhysMapHost_t >(_vec_map);
             std::get<0>(_borders) = *(_pm.getx().begin());
             std::get<1>(_borders) = *(_pm.getx().rbegin());
             std::get<2>(_borders) = *(_pm.gety().begin());
@@ -476,43 +489,35 @@ namespace maps{
             _multimap = false;
             std::cout << MIDDLEROW << std::endl;
             INFO_LINE("Volume boundaries: [ "<< std::get<0>(_borders) << " < x < "<< std::get<1>(_borders) << " ] micron, "<<"[ "<< std::get<2>(_borders) << " < y < "<< std::get<3>(_borders) <<" ] micron, "<<"[ "<< std::get<4>(_borders) << " < z < "<< std::get<5>(_borders) <<" ] micron")
+            _isloaded = true;
         }
 
         void load(std::string efield_path, std::string emob_path, std::string hmob_path, std::string wfield_path){
             INFO_LINE("Loading separate maps")
 
-            
+            //get efield
             _efieldm.init(efield_path); //E field is a vector
             _efieldm.prepare();
-            _efieldm.fillspacepoints<VecDev_t<double> >(_xvec_ef,_yvec_ef,_zvec_ef);
-            _efieldm.fillvector<VecDev_t<double> >(efieldvecx,efieldvecy,efieldvecz);
-            
-//             for(auto cf:configlist){
-//                 if(std::get<bool>(cf.get("plot"))){
-//                     VecHost_t<double> efieldvecx_host, efieldvecy_host, efieldvecz_host;
-//                     efieldm.fillvector<VecHost_t<double> >(efieldvecx_host,efieldvecy_host,efieldvecz_host);
-//                     h3=analysis::getHist3DDraw( efieldvecx_host, efieldvecy_host, efieldvecz_host, efieldm.getx(),efieldm.gety(),efieldm.getz());
-//                     break;
-//                 }
-//             }
+            _efieldm.fillspacepoints(_xvec_ef,_yvec_ef,_zvec_ef);
+            _efieldm.fillvector(_vec_ef);
 
 
             _wfieldm.init(wfield_path); // Weighting field is a vector
             _wfieldm.prepare();
-            _wfieldm.fillspacepoints<VecDev_t<double> >(_xvec_wf,_yvec_wf,_zvec_wf);
-            _wfieldm.fillvector<VecDev_t<double> >(wfieldvecx,wfieldvecy,wfieldvecz);
+            _wfieldm.fillspacepoints(_xvec_wf,_yvec_wf,_zvec_wf);
+            _wfieldm.fillvector(_vec_wf);
 
 
             _emobm.init(emob_path); // Electron mobility is a scalar
             _emobm.prepare();
-            _emobm.fillspacepoints<VecDev_t<double> >(_xvec_emob,_yvec_emob,_zvec_emob);
-            _emobm.fillscalar<VecDev_t<double> >(emobvec);
+            _emobm.fillspacepoints(_xvec_emob,_yvec_emob,_zvec_emob);
+            _emobm.fillscalar(emobvec);
 
 
             _hmobm.init(hmob_path); // Hole mobility is a scalar
             _hmobm.prepare();
-            _hmobm.fillspacepoints<VecDev_t<double> >(_xvec_hmob,_yvec_hmob,_zvec_hmob);
-            _hmobm.fillscalar<VecDev_t<double> >(hmobvec);
+            _hmobm.fillspacepoints(_xvec_hmob,_yvec_hmob,_zvec_hmob);
+            _hmobm.fillscalar(hmobvec);
 
 
             std::get<0>(_borders) = std::min(*(_efieldm.getx().begin()),std::min(*(_wfieldm.getx().begin()),std::min(*(_emobm.getx().begin()),*(_hmobm.getx().begin()))));
@@ -526,10 +531,13 @@ namespace maps{
 
             std::cout << MIDDLEROW << std::endl;
             INFO_LINE("Volume boundaries: [ "<< std::get<0>(_borders) << " < x < "<< std::get<1>(_borders) << " ] micron, "<<"[ "<< std::get<2>(_borders) << " < y < "<< std::get<3>(_borders) <<" ] micron, "<<"[ "<< std::get<4>(_borders) << " < z < "<< std::get<5>(_borders) <<" ] micron")
+            
+            _isloaded = true;
         }
 
         //is multimap or not?
         bool ismulti(){return _multimap;}
+        bool isloaded(){return _isloaded;}
     };
 
     //prepare maps

include/simulation/Evolve.h

@@ -227,12 +227,8 @@ namespace evolve{
 
             //calculate diffusion velocity
             double s = sqrt(2*(_fDiff*mobility)/_fStep)*hydra::get<_tc_gauss_x>(p);
-            hydra::Vector3R ref(1.,0.,0.);
 
-            double theta = hydra::get<_tc_angle_2>(p);
-            double phi = hydra::get<_tc_angle_1>(p);
-            
-            ref.set(s*sin(theta)*cos(phi),s*sin(theta)*sin(phi),s*cos(theta));
+            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)));
 
             //set diffusion velocity to be used in evolution step
             hydra::get<_tc_gauss_x>(p) = ref.get(0);
@@ -363,12 +359,8 @@ namespace evolve{
 
             //calculate diffusion velocity
             double s = sqrt(2*(_fDiff*mobility)/_fStep)*hydra::get<_tc_gauss_x>(p);
-            hydra::Vector3R ref(1.,0.,0.);
-            
-            double theta = hydra::get<_tc_angle_2>(p);
-            double phi = hydra::get<_tc_angle_1>(p);
 
-            ref.set(s*sin(theta)*cos(phi),s*sin(theta)*sin(phi),s*cos(theta));
+            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)));
 
             //set diffusion velocity to be used in evolution step
             hydra::get<_tc_gauss_x>(p) = ref.get(0);

include/simulation/HydraSim.h

@@ -45,19 +45,21 @@ void signal_simulation(cfg::Setting const& root, int icfg=-1){
 
 //     size_t nl=countlines(physmap_path);
 
+    
+    
     //3d hist of efield
     TH3D* h3=NULL;
 
     VecDev_t<double> _xvec, _yvec, _zvec;
-    hydra::multiarray<double,8,hydra::host::sys_t> _vec_map;
+    hydra::multiarray<double,8,hydra::device::sys_t> _vec_map;
 
     //efield vecs
     VecDev_t<double> _xvec_ef, _yvec_ef, _zvec_ef;
-    hydra::multiarray<double,3,hydra::host::sys_t> _vec_ef;
+    hydra::multiarray<double,3,hydra::device::sys_t> _vec_ef;
 
     //wfield vecs
     VecDev_t<double> _xvec_wf, _yvec_wf, _zvec_wf;
-    hydra::multiarray<double,3,hydra::host::sys_t> _vec_wf;
+    hydra::multiarray<double,3,hydra::device::sys_t> _vec_wf;
 
 
     //emob vecs
@@ -83,7 +85,7 @@ void signal_simulation(cfg::Setting const& root, int icfg=-1){
         pm.prepare();
 
         pm.fillspacepoints<VecDev_t<double> >(_xvec,_yvec,_zvec);
-        pm.fillallquantities<hydra::multiarray<double,8,hydra::host::sys_t> >(_vec_map);
+        pm.fillallquantities<hydra::multiarray<double,8,hydra::device::sys_t> >(_vec_map);
 
         _xmin = *(pm.getx().begin());
         _xmax = *(pm.getx().rbegin());
@@ -111,13 +113,13 @@ void signal_simulation(cfg::Setting const& root, int icfg=-1){
             maps::physmap<double,6> efieldm(efieldmap_path); //E field is a vector
             efieldm.prepare();
 
-            efieldm.fillspacepoints<VecDev_t<double> >(_xvec_ef,_yvec_ef,_zvec_ef);
-            efieldm.fillvector< hydra::multiarray<double,3,hydra::host::sys_t> >(_vec_ef);
+            efieldm.fillspacepoints(_xvec_ef,_yvec_ef,_zvec_ef);
+            efieldm.fillvector(_vec_ef);
 
             for(auto cf:configlist){
                 if(std::get<bool>(cf.get("plot"))){
                     VecHost_t<double> efieldvecx_host, efieldvecy_host, efieldvecz_host;
-                    efieldm.fillvector<VecHost_t<double> >(efieldvecx_host,efieldvecy_host,efieldvecz_host);
+                    efieldm.fillvector(efieldvecx_host,efieldvecy_host,efieldvecz_host);
                     h3=analysis::getHist3DDraw( efieldvecx_host, efieldvecy_host, efieldvecz_host, efieldm.getx(),efieldm.gety(),efieldm.getz());
                     break;
                 }
@@ -127,20 +129,20 @@ void signal_simulation(cfg::Setting const& root, int icfg=-1){
             std::string wfieldmap_path=root["PhysicsMaps"]["wfield"];
             maps::physmap<double,6> wfieldm(wfieldmap_path); // Weighting field is a vector
             wfieldm.prepare();
-            wfieldm.fillspacepoints<VecDev_t<double> >(_xvec_wf,_yvec_wf,_zvec_wf);
-            wfieldm.fillvector< hydra::multiarray<double,3,hydra::host::sys_t> >(_vec_wf);
+            wfieldm.fillspacepoints(_xvec_wf,_yvec_wf,_zvec_wf);
+            wfieldm.fillvector(_vec_wf);
 
             std::string emobmap_path=root["PhysicsMaps"]["emob"];
             maps::physmap<double,4> emobm(emobmap_path); // Electron mobility is a scalar
             emobm.prepare();
-            emobm.fillspacepoints<VecDev_t<double> >(_xvec_emob,_yvec_emob,_zvec_emob);
-            emobm.fillscalar<VecDev_t<double> >(emobvec);
+            emobm.fillspacepoints(_xvec_emob,_yvec_emob,_zvec_emob);
+            emobm.fillscalar(emobvec);
 
             std::string hmobmap_path=root["PhysicsMaps"]["hmob"];
             maps::physmap<double,4> hmobm(hmobmap_path); // Hole mobility is a scalar
             hmobm.prepare();
-            hmobm.fillspacepoints<VecDev_t<double> >(_xvec_hmob,_yvec_hmob,_zvec_hmob);
-            hmobm.fillscalar<VecDev_t<double> >(hmobvec);
+            hmobm.fillspacepoints(_xvec_hmob,_yvec_hmob,_zvec_hmob);
+            hmobm.fillscalar(hmobvec);
 
 
             _xmin = std::max(*(efieldm.getx().begin()),std::max(*(wfieldm.getx().begin()),std::max(*(emobm.getx().begin()),*(hmobm.getx().begin()))));