diff --git a/README.md b/README.md index 5a47bed..38305da 100644 --- a/README.md +++ b/README.md @@ -1,14 +1,14 @@ # TCode ## What is it? -TCode is a C++14 compliant application to simulate the response of solid state sensors in massively parallel platforms on Linux systems. TCode is implemented on top of [Hydra](https://github.com/MultithreadCorner/Hydra) and as such, it can run on OpenMP, CUDA and TBB compatible devices. +TCode is a C++14 compliant application to simulate the response of solid state sensors in massively parallel platforms on Linux systems. TCode is implemented on top of [Hydra](https://github.com/MultithreadCorner/Hydra) and as such, it can run in parallel on OpenMP, CUDA and TBB compatible devices or sequentialy in single threaded enviroments. TCode is still in its alpha version and the repository it is taking shape, for the moment we focused on getting raw performance and will make more stable and user friendly release in the near future. ## How it works -TCode uses external 3D maps of electric fields, carrier mobilities and weighting field and energy deposit to simulate the response in current of solid state sensors. The motion of the individual carriers produced in the initial deposit is determined using a 4th other Runge-Kutta using the electric field and the mobilities and assuming that the carriers always move at drift velocity. At each time interval the current induced in the electrod is calculated from the carriers velocity using the corresponding weighting field, according to the [Shockley](https://aip.scitation.org/doi/10.1063/1.1710367)-[Ramo](https://ieeexplore.ieee.org/document/1686997) theorem. The output is stored to ROOT files, several level of output detail are available, from the simple current vs time plot to the complete information about the position of the carriers at each time step. +TCode uses external 3D maps of electric fields, carrier mobilities and weighting field and energy deposit to simulate the response in current of solid state sensors. The motion of the individual carriers produced in the initial deposit is determined using a 4th oder Runge-Kutta method **[add reference]** using the electric field and the mobilities and assuming that the carriers always move at drift velocity. At regular times intervals, the current induced in the electron is calculated from the carriers velocity using the corresponding weighting field, according to the [Shockley](https://aip.scitation.org/doi/10.1063/1.1710367)-[Ramo](https://ieeexplore.ieee.org/document/1686997) theorem. The output is stored to ROOT files **(consider supporting flat text files as well )**, several level of output detail are available, from the simple current vs time plot to the complete information about the position of the carriers at each time step. ## Dependencies -TCode depends on [ROOT >= v.6.14](https://github.com/root-project/root), [libconfig >= v1.5](https://hyperrealm.github.io/libconfig/), [TCLAP >= v1.2.1](http://tclap.sourceforge.net/) and optionally [CUDA >= 10.0](https://developer.nvidia.com/cuda-toolkit) (needed for nVidia GPUs). +TCode depends on [ROOT >= v.6.14](https://github.com/root-project/root), [libconfig >= v1.5](https://hyperrealm.github.io/libconfig/), [TCLAP >= v1.2.1](http://tclap.sourceforge.net/) and optionally [CUDA >= 10.0](https://developer.nvidia.com/cuda-toolkit) (needed for nVidia GPUs) and TBB. **Perche cuda 10+ e non 8+?** ## Manual IN PREPARATION diff --git a/examples/example_config.cfg b/examples/example_config.cfg index 9ebdf06..fb7021a 100644 --- a/examples/example_config.cfg +++ b/examples/example_config.cfg @@ -46,7 +46,7 @@ InputData: #number of time steps in the simulation, integer. #NEEDED #it can be an array (e.g. [400,600]). In this case simulations will be performed with all the values - steps = 100; + steps = 600; #time step in seconds, double. NEEDED #it can be an array (e.g. [1e-11,1e-13,1e-14]). In this case simulations will be performed with all the values @@ -108,7 +108,7 @@ InputData: #it can be an array particles = 20000; #it cannot be -1 in this case since no input file is provided - T=0.; + T=300.; Xshift = 2.; Yshift = 2.; diff --git a/include/Evolve.h b/include/Evolve.h index 77e172c..3c5eda1 100644 --- a/include/Evolve.h +++ b/include/Evolve.h @@ -1,6 +1,6 @@ /*---------------------------------------------------------------------------- * - * Copyright (C) 2018 Andrea Contu e Angelo Loi + * Copyright (C) 2018-2019 Andrea Contu e Angelo Loi * * This file is part of TCode software. * @@ -22,9 +22,12 @@ * Evolve.h * * Created on: 12/11/2018 - * Author: Andrea Contu + * Author: Andrea Contu */ +#ifndef __EVOLVE_H__ +#define __EVOLVE_H__ + namespace evolve{ //functions to calculates trilinear interpolation @@ -44,11 +47,19 @@ namespace evolve{ } + /* + review A.A. + why not using std::binary_search (https://en.cppreference.com/w/cpp/algorithm/binary_search) ? + if V is too large (O(6) or higher) and it becomes too slow, we can wrapp a vectorized version + in hydra... + */ + //function to get index in one dimension (binary search) __hydra_dual__ inline size_t getindex(double v[], size_t n, double f, size_t first=0){ if(f<=v[first]){ return first+1; + } if(f>=v[n-1]){ return n-1; @@ -102,297 +113,28 @@ namespace evolve{ ind8[6] = ind[1]-1 +(ind[0])*sy+(ind[2])*fNPXY; ind8[7] = ind[1] + (ind[0])*sy+(ind[2])*fNPXY; } + /* + Review A.A. - //functor to fill carrier drift velocities - struct SetVeldiff - { - - SetVeldiff()=delete; // avoid creating objects with undefined data - - // No __hydra__dual__ here. Only allow objects to be constructed in host side. - SetVeldiff(size_t n, double timestep, double diffusion, double *in_x, size_t dim_x,double *in_y, size_t dim_y,double *in_z, size_t dim_z, double *in_em, double *in_hm): - fN(n), - fStep(timestep), - fDiff(diffusion), - fNPX(dim_x), - fNPY(dim_y), - fNPZ(dim_z), - xvec(in_x), - yvec(in_y), - zvec(in_z), - vecem(in_em), - vechm(in_hm) - {} - - // No __hydra__dual__ here. Only allow objects to be constructed in host side. - SetVeldiff(size_t n, double timestep, double diffusion, VecDev_t &in_x, VecDev_t &in_y, VecDev_t &in_z, VecDev_t &in_em, VecDev_t &in_hm): - fN(n), - fStep(timestep), - fDiff(diffusion), - fNPX(in_x.size()), - fNPY(in_y.size()), - fNPZ(in_z.size()), - xvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_x.data())), - yvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_y.data())), - zvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_z.data())), - vecem(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_em.data())), - vechm(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_hm.data())) - {} - - // Copy constructor needs be __hydra_dual__ to be copyable - // on device and host sides, in the time of thread distribution. - // hydra, as well as thrust and std pass arguments by value. - // Obs.: RooFit passes by reference... (crazy design choice!) - __hydra_dual__ - SetVeldiff( SetVeldiff const& other): - fN(other.fN), - fStep(other.fStep), - fDiff(other.fDiff), - fNPX(other.fNPX), - fNPY(other.fNPY), - fNPZ(other.fNPZ), - xvec(other.xvec), - yvec(other.yvec), - zvec(other.zvec), - vecem(other.vecem), - vechm(other.vechm) - { } - - - __hydra_dual__ - SetVeldiff operator=( SetVeldiff const& other){ - if(this == &other) return *this; - - fN=other.fN; - fStep=other.fStep; - fDiff=other.fDiff; - fNPX=other.fNPX; - fNPY=other.fNPY; - fNPZ=other.fNPZ; - xvec=other.xvec; - yvec=other.yvec; - zvec=other.zvec; - vecem=other.vecem; - vechm=other.vechm; - - return *this; - } - - //function call operator needs to be callable in host and device sides - // and be constant - template - __hydra_dual__ - void operator()(Particle p){ - - //check if particle went outside the volume - if(hydra::get<_tc_isin>(p)<0) return; - - //get particle charge and initial position - double charge=hydra::get<_tc_charge>(p); - double mobility = 0; - double x=hydra::get<_tc_x>(p); - double y=hydra::get<_tc_y>(p); - double z=hydra::get<_tc_z>(p); - - double vvec[3]; - size_t ind8[8]; - getallindices(vvec, ind8, xvec, fNPX, yvec, fNPY, zvec, fNPZ, x, y, z); - - if(charge<0.) mobility=interpol(vvec, vecem, ind8); - else mobility=interpol(vvec, vechm, ind8); - - //calculate diffusion - double D = ::sqrt(2*(fDiff*mobility)*fStep); //D=radice(2*kb*T*mu*dt/e) - double sx = D*hydra::get<_tc_gauss_x>(p); - double sy = D*hydra::get<_tc_gauss_y>(p); - double sz = D*hydra::get<_tc_gauss_z>(p); - - //spostamento in x - double alpha=hydra::get<_tc_angle_1>(p); - double beta=hydra::get<_tc_angle_2>(p); - double ax = cos(alpha)*cos(beta)*sx - sin(alpha)*sy - cos(alpha)*sin(beta)*sz; - //spostamento in y - double ay = sin(alpha)*cos(beta)*sx + cos(alpha)*sy - sin(alpha)*sin(beta)*sz; - //spostamento in z - double az = sin(beta)*sz + cos(beta)*sz; - - hydra::get<_tc_gauss_x>(p) = ax/fStep; - hydra::get<_tc_gauss_y>(p) = ay/fStep; - hydra::get<_tc_gauss_z>(p) = az/fStep; - } - - - size_t fN; - double fStep; - double fDiff; - size_t fNPX = 0; - size_t fNPY = 0; - size_t fNPZ = 0; - double *vecem = nullptr; - double *vechm = nullptr; - double *xvec = nullptr; - double *yvec = nullptr; - double *zvec = nullptr; - }; - - - //functor to fill carrier drift velocities - struct SetVeldiff_multi - { - - SetVeldiff_multi()=delete; // avoid creating objects with undefined data - - // No __hydra__dual__ here. Only allow objects to be constructed in host side. - SetVeldiff_multi(size_t n, double timestep, double diffusion, VecDev_t &ine_x, VecDev_t &ine_y, VecDev_t &ine_z, VecDev_t &inh_x, VecDev_t &inh_y, VecDev_t &inh_z, VecDev_t &in_em, VecDev_t &in_hm): - fN(n), - fStep(timestep), - fDiff(diffusion), - fNPXe(ine_x.size()), - fNPYe(ine_y.size()), - fNPZe(ine_z.size()), - xvece(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(ine_x.data())), - yvece(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(ine_y.data())), - zvece(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(ine_z.data())), - fNPXh(inh_x.size()), - fNPYh(inh_y.size()), - fNPZh(inh_z.size()), - xvech(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inh_x.data())), - yvech(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inh_y.data())), - zvech(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inh_z.data())), - vecem(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_em.data())), - vechm(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_hm.data())) - {} - - // Copy constructor needs be __hydra_dual__ to be copyable - // on device and host sides, in the time of thread distribution. - // hydra, as well as thrust and std pass arguments by value. - // Obs.: RooFit passes by reference... (crazy design choice!) - __hydra_dual__ - SetVeldiff_multi( SetVeldiff_multi const& other): - fN(other.fN), - fStep(other.fStep), - fDiff(other.fDiff), - fNPXe(other.fNPXe), - fNPYe(other.fNPYe), - fNPZe(other.fNPZe), - xvece(other.xvece), - yvece(other.yvece), - zvece(other.zvece), - fNPXh(other.fNPXh), - fNPYh(other.fNPYh), - fNPZh(other.fNPZh), - xvech(other.xvech), - yvech(other.yvech), - zvech(other.zvech), - vecem(other.vecem), - vechm(other.vechm) - { } - - - __hydra_dual__ - SetVeldiff_multi operator=( SetVeldiff_multi const& other){ - if(this == &other) return *this; - - fN=other.fN; - fStep=other.fStep; - fDiff=other.fDiff; - fNPXe=other.fNPXe; - fNPYe=other.fNPYe; - fNPZe=other.fNPZe; - xvece=other.xvece; - yvece=other.yvece; - zvece=other.zvece; - fNPXh=other.fNPXh; - fNPYh=other.fNPYh; - fNPZh=other.fNPZh; - xvech=other.xvech; - yvech=other.yvech; - zvech=other.zvech; - vecem=other.vecem; - vechm=other.vechm; - - return *this; - } - - //function call operator needs to be callable in host and device sides - // and be constant - template - __hydra_dual__ - void operator()(Particle p){ - - //check if particle went outside the volume - if(hydra::get<_tc_isin>(p)<0) return; - - //get particle charge and initial position - double charge=hydra::get<_tc_charge>(p); - double mobility = 0; - double x=hydra::get<_tc_x>(p); - double y=hydra::get<_tc_y>(p); - double z=hydra::get<_tc_z>(p); - - double vvec[3]; - size_t ind8[8]; - - if(charge<0.){ - getallindices(vvec, ind8, xvece, fNPXe, yvece, fNPYe, zvece, fNPZe, x, y, z); - mobility=interpol(vvec, vecem, ind8); - } - else{ - getallindices(vvec, ind8, xvech, fNPXh, yvech, fNPYh, zvech, fNPZh, x, y, z); - mobility=interpol(vvec, vechm, ind8); - } - - //calculate diffusion - double D = ::sqrt(2*(fDiff*mobility)*fStep); //D=radice(2*kb*T*mu*dt/e) - double sx = D*hydra::get<_tc_gauss_x>(p); - double sy = D*hydra::get<_tc_gauss_y>(p); - double sz = D*hydra::get<_tc_gauss_z>(p); - - //spostamento in x - double alpha=hydra::get<_tc_angle_1>(p); - double beta=hydra::get<_tc_angle_2>(p); - double ax = cos(alpha)*cos(beta)*sx - sin(alpha)*sy - cos(alpha)*sin(beta)*sz; - //spostamento in y - double ay = sin(alpha)*cos(beta)*sx + cos(alpha)*sy - sin(alpha)*sin(beta)*sz; - //spostamento in z - double az = sin(beta)*sz + cos(beta)*sz; - - hydra::get<_tc_gauss_x>(p) = ax/fStep; - hydra::get<_tc_gauss_y>(p) = ay/fStep; - hydra::get<_tc_gauss_z>(p) = az/fStep; - } - - - size_t fN; - double fStep; - double fDiff; - size_t fNPXe = 0; - size_t fNPYe = 0; - size_t fNPZe = 0; - size_t fNPXh = 0; - size_t fNPYh = 0; - size_t fNPZh = 0; - double *vecem = nullptr; - double *vechm = nullptr; - double *xvece = nullptr; - double *yvece = nullptr; - double *zvece = nullptr; - double *xvech = nullptr; - double *yvech = nullptr; - double *zvech = nullptr; - }; - + ApplyRamo is receiving and using pointers without testing the validity. + */ //calculate currents struct ApplyRamo { ApplyRamo()=delete; // avoid creating objects with undefined data + /* + Review A.A. + + ApplyRamo ctor is receiving and using pointers without testing the validity. + */ // No __hydra__dual__ here. Only allow objects to be constructed in host side. - ApplyRamo(size_t n, double timestep, double *in_x, size_t dim_x,double *in_y, size_t dim_y,double *in_z, size_t dim_z, double *in_xef,double *in_yef, double *in_zef, double *in_em, double *in_hm, double *in_xwf,double *in_ywf, double *in_zwf): + ApplyRamo(size_t n, double timestep, double diffusion, double *in_x, size_t dim_x,double *in_y, size_t dim_y,double *in_z, size_t dim_z, double *in_xef,double *in_yef, double *in_zef, double *in_em, double *in_hm, double *in_xwf,double *in_ywf, double *in_zwf): fN(n), fStep(timestep), + fDiff(diffusion), fNPX(dim_x), fNPY(dim_y), fNPZ(dim_z), @@ -409,10 +151,12 @@ namespace evolve{ zvecwf(in_zwf) {} + // // No __hydra__dual__ here. Only allow objects to be constructed in host side. - ApplyRamo(size_t n, double timestep, VecDev_t &in_x, VecDev_t &in_y, VecDev_t &in_z, VecDev_t &in_xef, VecDev_t &in_yef, VecDev_t &in_zef, VecDev_t &in_em, VecDev_t &in_hm, VecDev_t &in_xwf, VecDev_t &in_ywf, VecDev_t &in_zwf): + ApplyRamo(size_t n, double timestep, double diffusion, VecDev_t &in_x, VecDev_t &in_y, VecDev_t &in_z, VecDev_t &in_xef, VecDev_t &in_yef, VecDev_t &in_zef, VecDev_t &in_em, VecDev_t &in_hm, VecDev_t &in_xwf, VecDev_t &in_ywf, VecDev_t &in_zwf): fN(n), fStep(timestep), + fDiff(diffusion), fNPX(in_x.size()), fNPY(in_y.size()), fNPZ(in_z.size()), @@ -437,6 +181,7 @@ namespace evolve{ ApplyRamo( ApplyRamo const& other): fN(other.fN), fStep(other.fStep), + fDiff(other.fDiff), fNPX(other.fNPX), fNPY(other.fNPY), fNPZ(other.fNPZ), @@ -459,6 +204,7 @@ namespace evolve{ fN=other.fN; fStep=other.fStep; + fDiff=other.fDiff; fNPX=other.fNPX; fNPY=other.fNPY; fNPZ=other.fNPZ; @@ -477,6 +223,7 @@ namespace evolve{ return *this; } + // FIXME: is not constant ! //function call operator needs to be callable in host and device sides // and be constant template @@ -518,9 +265,31 @@ namespace evolve{ double k1y=charge*mobility*ey; double k1z=charge*mobility*ez; - double Vdriftx=k1x+hydra::get<_tc_gauss_x>(p); - double Vdrifty=k1y+hydra::get<_tc_gauss_y>(p); - double Vdriftz=k1z+hydra::get<_tc_gauss_z>(p); + //calculate diffusion velocity + double s = ::sqrt(2*(fDiff*mobility)/fStep)*hydra::get<_tc_gauss_x>(p); + + hydra::Vector3R direction(ex,ey,ez); + hydra::Vector3R ref(1.,0.,0.); + double theta = hydra::get<_tc_angle_2>(p); + double phi = hydra::get<_tc_angle_1>(p); + if(direction.d3mag()==0){ + ref.set(s*sin(theta)*cos(phi),s*sin(theta)*sin(phi),s*cos(theta)); + } + else{ + + direction = direction.unit(); + ref = (direction.cross(ref)).unit()*s; + //rotate ref + ref = ref*cos(phi) + (direction.cross(ref))*sin(phi)+direction*(direction.dot(ref))*(1-cos(phi)); + } + + hydra::get<_tc_gauss_x>(p) = ref.get(0); + hydra::get<_tc_gauss_y>(p) = ref.get(1); + hydra::get<_tc_gauss_z>(p) = ref.get(2); + + double Vdriftx=k1x+ref.get(0); + double Vdrifty=k1y+ref.get(1); + double Vdriftz=k1z+ref.get(2); // Vdrift=sqrt((Vdriftx*Vdriftx+Vdrifty*Vdrifty+Vdriftz*Vdriftz)); @@ -539,9 +308,11 @@ namespace evolve{ } - + //FIXME: pointers initialized here and in the class ctors. Consider to enable the default ctor. + //anyway you never test the pointers for null. size_t fN; double fStep; + double fDiff; size_t fNPX = 0; size_t fNPY = 0; size_t fNPZ = 0; @@ -558,17 +329,20 @@ namespace evolve{ double *zvec = nullptr; }; - + //fixme: why not using the iterators instead of raw pointer? + // your only ctor is taking vectors ... //calculate currents struct ApplyRamo_multi { + //fixme: Actually you are anyway initliazing the struct with null-pointers ApplyRamo_multi()=delete; // avoid creating objects with undefined data // No __hydra__dual__ here. Only allow objects to be constructed in host side. - ApplyRamo_multi(size_t n, double timestep, VecDev_t &inef_x, VecDev_t &inef_y, VecDev_t &inef_z, VecDev_t &in_xef, VecDev_t &in_yef, VecDev_t &in_zef, VecDev_t &inem_x, VecDev_t &inem_y, VecDev_t &inem_z, VecDev_t &in_em, VecDev_t &inhm_x, VecDev_t &inhm_y, VecDev_t &inhm_z, VecDev_t &in_hm, VecDev_t &inwf_x, VecDev_t &inwf_y, VecDev_t &inwf_z, VecDev_t &in_xwf, VecDev_t &in_ywf, VecDev_t &in_zwf): + ApplyRamo_multi(size_t n, double timestep, double diffusion, VecDev_t &inef_x, VecDev_t &inef_y, VecDev_t &inef_z, VecDev_t &in_xef, VecDev_t &in_yef, VecDev_t &in_zef, VecDev_t &inem_x, VecDev_t &inem_y, VecDev_t &inem_z, VecDev_t &in_em, VecDev_t &inhm_x, VecDev_t &inhm_y, VecDev_t &inhm_z, VecDev_t &in_hm, VecDev_t &inwf_x, VecDev_t &inwf_y, VecDev_t &inwf_z, VecDev_t &in_xwf, VecDev_t &in_ywf, VecDev_t &in_zwf): fN(n), fStep(timestep), + fDiff(diffusion), fNPXef(inef_x.size()), fNPYef(inef_y.size()), fNPZef(inef_z.size()), @@ -611,6 +385,7 @@ namespace evolve{ ApplyRamo_multi( ApplyRamo_multi const& other): fN(other.fN), fStep(other.fStep), + fDiff(other.fDiff), fNPXef(other.fNPXef), fNPYef(other.fNPYef), fNPZef(other.fNPZef), @@ -651,6 +426,7 @@ namespace evolve{ fN=other.fN; fStep=other.fStep; + fDiff=other.fDiff; fNPXef=other.fNPXef; fNPYef=other.fNPYef; fNPZef=other.fNPZef; @@ -734,10 +510,31 @@ namespace evolve{ double k1x=charge*mobility*ex; double k1y=charge*mobility*ey; double k1z=charge*mobility*ez; + //calculate diffusion velocity + double s = sqrt(2*(fDiff*mobility)/fStep)*hydra::get<_tc_gauss_x>(p); + hydra::Vector3R direction(ex,ey,ez); + hydra::Vector3R ref(1.,0.,0.); + + double theta = hydra::get<_tc_angle_2>(p); + double phi = hydra::get<_tc_angle_1>(p); + if(direction.d3mag()==0){ + ref.set(s*sin(theta)*cos(phi),s*sin(theta)*sin(phi),s*cos(theta)); + } + else{ + + direction = direction.unit(); + ref = (direction.cross(ref)).unit()*s; + //rotate ref + ref = ref*cos(phi) + (direction.cross(ref))*sin(phi)+direction*(direction.dot(ref))*(1-cos(phi)); + } - double Vdriftx=k1x+hydra::get<_tc_gauss_x>(p); - double Vdrifty=k1y+hydra::get<_tc_gauss_y>(p); - double Vdriftz=k1z+hydra::get<_tc_gauss_z>(p); + hydra::get<_tc_gauss_x>(p) = ref.get(0); + hydra::get<_tc_gauss_y>(p) = ref.get(1); + hydra::get<_tc_gauss_z>(p) = ref.get(2); + + double Vdriftx=k1x+ref.get(0); + double Vdrifty=k1y+ref.get(1); + double Vdriftz=k1z+ref.get(2); // Vdrift=sqrt((Vdriftx*Vdriftx+Vdrifty*Vdrifty+Vdriftz*Vdriftz)); @@ -758,6 +555,7 @@ namespace evolve{ size_t fN; double fStep; + double fDiff; size_t fNPXef = 0; size_t fNPYef = 0; size_t fNPZef = 0; @@ -792,9 +590,10 @@ namespace evolve{ double *zhmvec = nullptr; }; + //fixme: test pointer validity before to use! struct Evolve { - + //fixme: Actually you are anyway initliazing the struct with null-pointers Evolve()=delete; // avoid creating objects with undefined data // No __hydra__dual__ here. Only allow objects to be constructed in host side. @@ -815,7 +614,7 @@ namespace evolve{ {} // No __hydra__dual__ here. Only allow objects to be constructed in host side. - Evolve(size_t n, double timestep, VecDev_t &in_x, VecDev_t &in_y, VecDev_t &in_z, VecDev_t &in_xef, VecDev_t &in_yef, VecDev_t &in_zef, VecDev_t &in_em, VecDev_t &in_hm): + Evolve(size_t n, double timestep, VecDev_t &in_x, VecDev_t &in_y, VecDev_t &in_z, VecDev_t &in_xef, VecDev_t &in_yef, VecDev_t &in_zef, VecDev_t &in_em, VecDev_t &in_hm): fN(n), fStep(timestep), fNPX(in_x.size()), @@ -1023,7 +822,7 @@ namespace evolve{ Evolve_multi()=delete; // avoid creating objects with undefined data // No __hydra__dual__ here. Only allow objects to be constructed in host side. - Evolve_multi(size_t n, double timestep, VecDev_t &inef_x, VecDev_t &inef_y, VecDev_t &inef_z, VecDev_t &in_xef, VecDev_t &in_yef, VecDev_t &in_zef, VecDev_t &inem_x, VecDev_t &inem_y, VecDev_t &inem_z, VecDev_t &in_em, VecDev_t &inhm_x, VecDev_t &inhm_y, VecDev_t &inhm_z, VecDev_t &in_hm, double _xmin, double _xmax, double _ymin, double _ymax, double _zmin, double _zmax): + Evolve_multi(size_t n, double timestep, VecDev_t &inef_x, VecDev_t &inef_y, VecDev_t &inef_z, VecDev_t &in_xef, VecDev_t &in_yef, VecDev_t &in_zef, VecDev_t &inem_x, VecDev_t &inem_y, VecDev_t &inem_z, VecDev_t &in_em, VecDev_t &inhm_x, VecDev_t &inhm_y, VecDev_t &inhm_z, VecDev_t &in_hm, double _xmin, double _xmax, double _ymin, double _ymax, double _zmin, double _zmax): fN(n), fStep(timestep), xmin(_xmin), @@ -1325,3 +1124,4 @@ namespace evolve{ double *zhmvec = nullptr; }; } +#endif diff --git a/include/Evolve_dev.h b/include/Evolve_dev.h deleted file mode 100644 index cedb9d3..0000000 --- a/include/Evolve_dev.h +++ /dev/null @@ -1,1390 +0,0 @@ -/*---------------------------------------------------------------------------- - * - * Copyright (C) 2018 Andrea Contu e Angelo Loi - * - * This file is part of SimSens software. - * - * SimSens is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * SimSens is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with Hydra. If not, see . - * - *---------------------------------------------------------------------------*/ -/* - * Evolve.h - * - * Created on: 12/11/2018 - * Author: Andrea Contu - */ - -namespace evolve{ - - //functions to rotate around vector - __hydra_dual__ - inline void _getdiffvect(double vecdir[],double startvec[], double outvec[], double angle, double modulus){ - double cosphi = cos(angle); - double sinphi = sin(angle); - double onemcosphi = 1 - cosphi; - - //define rotation matric around a vector - double R11 = cosphi+vecdir[0]*vecdir[0]*onemcosphi; - double R21 = vecdir[0]*vecdir[1]*onemcosphi+vecdir[2]*sinphi; - double R31 = vecdir[0]*vecdir[2]*onemcosphi-vecdir[1]*sinphi; - - double R12 = vecdir[0]*vecdir[1]*onemcosphi-vecdir[2]*sinphi; - double R22 = cosphi+vecdir[1]*vecdir[1]*onemcosphi; - double R32 = vecdir[1]*vecdir[2]*onemcosphi+vecdir[0]*sinphi; - - double R13 = vecdir[0]*vecdir[2]*onemcosphi+vecdir[1]*sinphi; - double R23 = vecdir[1]*vecdir[2]*onemcosphi-vecdir[0]*sinphi; - double R33 = cosphi+vecdir[2]*vecdir[2]*onemcosphi; - - outvec[0] = (startvec[0]*R11+startvec[1]*R12+startvec[2]*R13)*modulus; - outvec[1] = (startvec[0]*R21+startvec[1]*R22+startvec[2]*R23)*modulus; - outvec[0] = (startvec[0]*R31+startvec[1]*R32+startvec[2]*R33)*modulus; - } - - __hydra_dual__ - inline void getdiffvect(double vecdir[], double outvec[], double angle, double modulus){ - double _vecdir[3]; - double modvecdir = sqrt(vecdir[0]*vecdir[0]+vecdir[1]*vecdir[1]+vecdir[2]*vecdir[2]); - _vecdir[0] = vecdir[0]/modvecdir; - _vecdir[1] = vecdir[1]/modvecdir; - _vecdir[2] = vecdir[2]/modvecdir; - double startvec[3]; - startvec[2] = 0.; - startvec[1] = 1.; - startvec[0] = -vecdir[1]/vecdir[0]*startvec[1]; - double startvecmod = sqrt(startvec[0]*startvec[0]+startvec[1]*startvec[1]+startvec[2]*startvec[2]); - startvec[0] = startvec[0]/startvecmod; - startvec[1] = startvec[1]/startvecmod; - startvec[2] = startvec[2]/startvecmod; - - _getdiffvect(_vecdir, startvec, outvec, angle, modulus); - } - - //functions to calculates trilinear interpolation - __hydra_dual__ - inline double cfun(double d, double t1, double t2){ - return t1*(1.-d)+t2*d; - } - - __hydra_dual__ - inline double val(double xd, double yd, double zd, double p0, double px, double py, double pz, double pxy, double pyz, double pxz, double pxyz){ - return cfun(zd,cfun(yd, cfun(xd,p0,px),cfun(xd,py,pxy)), cfun(yd,cfun(xd,pz,pxz), cfun(xd,pyz,pxyz))); - } - - __hydra_dual__ - inline double interpol(double vvec[], double v[], size_t ind8[8]){ - return val(vvec[0], vvec[1], vvec[2], v[ind8[0]], v[ind8[1]], v[ind8[2]], v[ind8[3]], v[ind8[4]], v[ind8[5]], v[ind8[6]], v[ind8[7]]); - - } - - //function to get index in one dimension (binary search) - __hydra_dual__ - inline size_t getindex(double v[], size_t n, double f, size_t first=0){ - if(f<=v[first]){ - return first+1; - } - if(f>=v[n-1]){ - return n-1; - } - size_t lo = first; - size_t hi = n - 1; - size_t mid=0; - while(lo0 && f>v[mid-1]) return mid; - hi=mid; - } - else{ - if(mid - __hydra_dual__ - void operator()(Particle p){ - - //check if particle went outside the volume - if(hydra::get<_tc_isin>(p)<0) return; - - //get particle charge and initial position - double charge=hydra::get<_tc_charge>(p); - double mobility = 0; - double x=hydra::get<_tc_x>(p); - double y=hydra::get<_tc_y>(p); - double z=hydra::get<_tc_z>(p); - - double vvec[3]; - size_t ind8[8]; - getallindices(vvec, ind8, xvec, fNPX, yvec, fNPY, zvec, fNPZ, x, y, z); - - if(charge<0.) mobility=interpol(vvec, vecem, ind8); - else mobility=interpol(vvec, vechm, ind8); - - //calculate diffusion - double D = ::sqrt(2*(fDiff*mobility)*fStep); //D=radice(2*kb*T*mu*dt/e) - double sx = D*hydra::get<_tc_gauss_x>(p); -// double sy = D*hydra::get<_tc_gauss_y>(p); -// double sz = D*hydra::get<_tc_gauss_z>(p); - - - //spostamento in x -// double alpha=hydra::get<_tc_angle_1>(p); -// double beta=hydra::get<_tc_angle_2>(p); -// double ax = cos(alpha)*cos(beta)*sx - sin(alpha)*sy - cos(alpha)*sin(beta)*sz; -// //spostamento in y -// double ay = sin(alpha)*cos(beta)*sx + cos(alpha)*sy - sin(alpha)*sin(beta)*sz; -// //spostamento in z -// double az = sin(beta)*sz + cos(beta)*sz; - - hydra::get<_tc_gauss_x>(p) = sx; -// hydra::get<_tc_gauss_y>(p) = ay/fStep; -// hydra::get<_tc_gauss_z>(p) = az/fStep; - } - - - size_t fN; - double fStep; - double fDiff; - size_t fNPX = 0; - size_t fNPY = 0; - size_t fNPZ = 0; - double *vecem = nullptr; - double *vechm = nullptr; - double *xvec = nullptr; - double *yvec = nullptr; - double *zvec = nullptr; - }; - - - //functor to fill carrier drift velocities - struct SetVeldiff_multi - { - - SetVeldiff_multi()=delete; // avoid creating objects with undefined data - - // No __hydra__dual__ here. Only allow objects to be constructed in host side. - SetVeldiff_multi(size_t n, double timestep, double diffusion, VecDev_t &ine_x, VecDev_t &ine_y, VecDev_t &ine_z, VecDev_t &inh_x, VecDev_t &inh_y, VecDev_t &inh_z, VecDev_t &in_em, VecDev_t &in_hm): - fN(n), - fStep(timestep), - fDiff(diffusion), - fNPXe(ine_x.size()), - fNPYe(ine_y.size()), - fNPZe(ine_z.size()), - xvece(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(ine_x.data())), - yvece(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(ine_y.data())), - zvece(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(ine_z.data())), - fNPXh(inh_x.size()), - fNPYh(inh_y.size()), - fNPZh(inh_z.size()), - xvech(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inh_x.data())), - yvech(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inh_y.data())), - zvech(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inh_z.data())), - vecem(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_em.data())), - vechm(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_hm.data())) - {} - - // Copy constructor needs be __hydra_dual__ to be copyable - // on device and host sides, in the time of thread distribution. - // hydra, as well as thrust and std pass arguments by value. - // Obs.: RooFit passes by reference... (crazy design choice!) - __hydra_dual__ - SetVeldiff_multi( SetVeldiff_multi const& other): - fN(other.fN), - fStep(other.fStep), - fDiff(other.fDiff), - fNPXe(other.fNPXe), - fNPYe(other.fNPYe), - fNPZe(other.fNPZe), - xvece(other.xvece), - yvece(other.yvece), - zvece(other.zvece), - fNPXh(other.fNPXh), - fNPYh(other.fNPYh), - fNPZh(other.fNPZh), - xvech(other.xvech), - yvech(other.yvech), - zvech(other.zvech), - vecem(other.vecem), - vechm(other.vechm) - { } - - - __hydra_dual__ - SetVeldiff_multi operator=( SetVeldiff_multi const& other){ - if(this == &other) return *this; - - fN=other.fN; - fStep=other.fStep; - fDiff=other.fDiff; - fNPXe=other.fNPXe; - fNPYe=other.fNPYe; - fNPZe=other.fNPZe; - xvece=other.xvece; - yvece=other.yvece; - zvece=other.zvece; - fNPXh=other.fNPXh; - fNPYh=other.fNPYh; - fNPZh=other.fNPZh; - xvech=other.xvech; - yvech=other.yvech; - zvech=other.zvech; - vecem=other.vecem; - vechm=other.vechm; - - return *this; - } - - //function call operator needs to be callable in host and device sides - // and be constant - template - __hydra_dual__ - void operator()(Particle p){ - - //check if particle went outside the volume - if(hydra::get<_tc_isin>(p)<0) return; - - //get particle charge and initial position - double charge=hydra::get<_tc_charge>(p); - double mobility = 0; - double x=hydra::get<_tc_x>(p); - double y=hydra::get<_tc_y>(p); - double z=hydra::get<_tc_z>(p); - - double vvec[3]; - size_t ind8[8]; - - if(charge<0.){ - getallindices(vvec, ind8, xvece, fNPXe, yvece, fNPYe, zvece, fNPZe, x, y, z); - mobility=interpol(vvec, vecem, ind8); - } - else{ - getallindices(vvec, ind8, xvech, fNPXh, yvech, fNPYh, zvech, fNPZh, x, y, z); - mobility=interpol(vvec, vechm, ind8); - } - - //calculate diffusion - double D = ::sqrt(2*(fDiff*mobility)*fStep); //D=radice(2*kb*T*mu*dt/e) - double sx = D*hydra::get<_tc_gauss_x>(p); -// double sy = D*hydra::get<_tc_gauss_y>(p); -// double sz = D*hydra::get<_tc_gauss_z>(p); - - //spostamento in x -// double alpha=hydra::get<_tc_angle_1>(p); -// double beta=hydra::get<_tc_angle_2>(p); -// double ax = cos(alpha)*cos(beta)*sx - sin(alpha)*sy - cos(alpha)*sin(beta)*sz; -// //spostamento in y -// double ay = sin(alpha)*cos(beta)*sx + cos(alpha)*sy - sin(alpha)*sin(beta)*sz; -// //spostamento in z -// double az = sin(beta)*sz + cos(beta)*sz; - - hydra::get<_tc_gauss_x>(p) = sx; -// hydra::get<_tc_gauss_y>(p) = ay/fStep; -// hydra::get<_tc_gauss_z>(p) = az/fStep; - } - - - size_t fN; - double fStep; - double fDiff; - size_t fNPXe = 0; - size_t fNPYe = 0; - size_t fNPZe = 0; - size_t fNPXh = 0; - size_t fNPYh = 0; - size_t fNPZh = 0; - double *vecem = nullptr; - double *vechm = nullptr; - double *xvece = nullptr; - double *yvece = nullptr; - double *zvece = nullptr; - double *xvech = nullptr; - double *yvech = nullptr; - double *zvech = nullptr; - }; - - - //calculate currents - struct ApplyRamo - { - - ApplyRamo()=delete; // avoid creating objects with undefined data - - // No __hydra__dual__ here. Only allow objects to be constructed in host side. - ApplyRamo(size_t n, double timestep, double *in_x, size_t dim_x,double *in_y, size_t dim_y,double *in_z, size_t dim_z, double *in_xef,double *in_yef, double *in_zef, double *in_em, double *in_hm, double *in_xwf,double *in_ywf, double *in_zwf): - fN(n), - fStep(timestep), - fNPX(dim_x), - fNPY(dim_y), - fNPZ(dim_z), - xvec(in_x), - yvec(in_y), - zvec(in_z), - xvecef(in_xef), - yvecef(in_yef), - zvecef(in_zef), - vecem(in_em), - vechm(in_hm), - xvecwf(in_xwf), - yvecwf(in_ywf), - zvecwf(in_zwf) - {} - - // No __hydra__dual__ here. Only allow objects to be constructed in host side. - ApplyRamo(size_t n, double timestep, VecDev_t &in_x, VecDev_t &in_y, VecDev_t &in_z, VecDev_t &in_xef, VecDev_t &in_yef, VecDev_t &in_zef, VecDev_t &in_em, VecDev_t &in_hm, VecDev_t &in_xwf, VecDev_t &in_ywf, VecDev_t &in_zwf): - fN(n), - fStep(timestep), - fNPX(in_x.size()), - fNPY(in_y.size()), - fNPZ(in_z.size()), - xvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_x.data())), - yvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_y.data())), - zvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_z.data())), - xvecef(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_xef.data())), - yvecef(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_yef.data())), - zvecef(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_zef.data())), - vecem(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_em.data())), - vechm(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_hm.data())), - xvecwf(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_xwf.data())), - yvecwf(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_ywf.data())), - zvecwf(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_zwf.data())) - {} - - // Copy constructor needs be __hydra_dual__ to be copyable - // on device and host sides, in the time of thread distribution. - // hydra, as well as thrust and std pass arguments by value. - // Obs.: RooFit passes by reference... (crazy design choice!) - __hydra_dual__ - ApplyRamo( ApplyRamo const& other): - fN(other.fN), - fStep(other.fStep), - fNPX(other.fNPX), - fNPY(other.fNPY), - fNPZ(other.fNPZ), - xvec(other.xvec), - yvec(other.yvec), - zvec(other.zvec), - xvecef(other.xvecef), - yvecef(other.yvecef), - zvecef(other.zvecef), - vecem(other.vecem), - vechm(other.vechm), - xvecwf(other.xvecwf), - yvecwf(other.yvecwf), - zvecwf(other.zvecwf) - { } - - __hydra_dual__ - ApplyRamo operator=( ApplyRamo const& other){ - if(this == &other) return *this; - - fN=other.fN; - fStep=other.fStep; - fNPX=other.fNPX; - fNPY=other.fNPY; - fNPZ=other.fNPZ; - xvec=other.xvec; - yvec=other.yvec; - zvec=other.zvec; - xvecef=other.xvecef; - yvecef=other.yvecef; - zvecef=other.zvecef; - vecem=other.vecem; - vechm=other.vechm; - xvecwf=other.xvecwf; - yvecwf=other.yvecwf; - zvecwf=other.zvecwf; - - return *this; - } - - //function call operator needs to be callable in host and device sides - // and be constant - template - __hydra_dual__ - void operator()(Particle p){ - - //check if particle went outside the volume - if(hydra::get<_tc_isin>(p)<0){ - hydra::get<_tc_curr>(p)=0.; - return; - } - - //get particle charge and initial position - double charge=hydra::get<0>(p); - double mobility=0; - double x=hydra::get<_tc_x>(p); - double y=hydra::get<_tc_y>(p); - double z=hydra::get<_tc_z>(p); - - double vvec[3]; - size_t ind8[8]; - - getallindices(vvec, ind8, xvec, fNPX, yvec, fNPY, zvec, fNPZ, x, y, z); - if(charge<0.) mobility=interpol(vvec, vecem, ind8); - else mobility=interpol(vvec, vechm, ind8); - - //get electric field - double ex = interpol(vvec, xvecef, ind8); - double ey = interpol(vvec, yvecef, ind8); - double ez = interpol(vvec, zvecef, ind8); - - //get weighting field - double wfx = interpol(vvec, xvecwf, ind8); - double wfy = interpol(vvec, yvecwf, ind8); - double wfz = interpol(vvec, zvecwf, ind8); - - //first step only - double k1x=charge*mobility*ex; - double k1y=charge*mobility*ey; - double k1z=charge*mobility*ez; - - //set diff velocity vector - double veldiff[3]={0,0,0}; - double velpoint[3]={k1x,k1y,k1z}; - getdiffvect(velpoint, veldiff, hydra::get<_tc_angle_1>(p), hydra::get<_tc_gauss_x>(p)); - - hydra::get<_tc_gauss_x>(p) = veldiff[0]; - hydra::get<_tc_gauss_y>(p) = veldiff[1]; - hydra::get<_tc_gauss_z>(p) = veldiff[2]; - - double Vdriftx=k1x+hydra::get<_tc_gauss_x>(p); - double Vdrifty=k1y+hydra::get<_tc_gauss_y>(p); - double Vdriftz=k1z+hydra::get<_tc_gauss_z>(p); - -// Vdrift=sqrt((Vdriftx*Vdriftx+Vdrifty*Vdrifty+Vdriftz*Vdriftz)); - -// MaxV=MAXVDRIFTE; -// if(charge>0) MaxV=MAXVDRIFTH; -// -// if(Vdrift>MaxV){ -// Vdriftx*=(MaxV/(Vdrift)); -// Vdrifty*=(MaxV/(Vdrift)); -// Vdriftz*=(MaxV/(Vdrift)); -// } - - //calculate current at this instant - hydra::get<5>(p)=charge*HCHARGE*((Vdriftx)*wfx+(Vdrifty)*wfy+(Vdriftz)*wfz); - - - } - - - size_t fN; - double fStep; - size_t fNPX = 0; - size_t fNPY = 0; - size_t fNPZ = 0; - double *xvecwf = nullptr; - double *yvecwf = nullptr; - double *zvecwf = nullptr; - double *xvecef = nullptr; - double *yvecef = nullptr; - double *zvecef = nullptr; - double *vecem = nullptr; - double *vechm = nullptr; - double *xvec = nullptr; - double *yvec = nullptr; - double *zvec = nullptr; - }; - - - //calculate currents - struct ApplyRamo_multi - { - - ApplyRamo_multi()=delete; // avoid creating objects with undefined data - - // No __hydra__dual__ here. Only allow objects to be constructed in host side. - ApplyRamo_multi(size_t n, double timestep, VecDev_t &inef_x, VecDev_t &inef_y, VecDev_t &inef_z, VecDev_t &in_xef, VecDev_t &in_yef, VecDev_t &in_zef, VecDev_t &inem_x, VecDev_t &inem_y, VecDev_t &inem_z, VecDev_t &in_em, VecDev_t &inhm_x, VecDev_t &inhm_y, VecDev_t &inhm_z, VecDev_t &in_hm, VecDev_t &inwf_x, VecDev_t &inwf_y, VecDev_t &inwf_z, VecDev_t &in_xwf, VecDev_t &in_ywf, VecDev_t &in_zwf): - fN(n), - fStep(timestep), - fNPXef(inef_x.size()), - fNPYef(inef_y.size()), - fNPZef(inef_z.size()), - xefvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inef_x.data())), - yefvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inef_y.data())), - zefvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inef_z.data())), - fNPXwf(inwf_x.size()), - fNPYwf(inwf_y.size()), - fNPZwf(inwf_z.size()), - xwfvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inwf_x.data())), - ywfvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inwf_y.data())), - zwfvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inwf_z.data())), - fNPXem(inem_x.size()), - fNPYem(inem_y.size()), - fNPZem(inem_z.size()), - xemvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inem_x.data())), - yemvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inem_y.data())), - zemvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inem_z.data())), - fNPXhm(inhm_x.size()), - fNPYhm(inhm_y.size()), - fNPZhm(inhm_z.size()), - xhmvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inhm_x.data())), - yhmvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inhm_y.data())), - zhmvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inhm_z.data())), - xvecef(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_xef.data())), - yvecef(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_yef.data())), - zvecef(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_zef.data())), - vecem(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_em.data())), - vechm(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_hm.data())), - xvecwf(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_xwf.data())), - yvecwf(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_ywf.data())), - zvecwf(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_zwf.data())) - {} - - // Copy constructor needs be __hydra_dual__ to be copyable - // on device and host sides, in the time of thread distribution. - // hydra, as well as thrust and std pass arguments by value. - // Obs.: RooFit passes by reference... (crazy design choice!) - __hydra_dual__ - ApplyRamo_multi( ApplyRamo_multi const& other): - fN(other.fN), - fStep(other.fStep), - fNPXef(other.fNPXef), - fNPYef(other.fNPYef), - fNPZef(other.fNPZef), - xefvec(other.xefvec), - yefvec(other.yefvec), - zefvec(other.zefvec), - fNPXwf(other.fNPXwf), - fNPYwf(other.fNPYwf), - fNPZwf(other.fNPZwf), - xwfvec(other.xwfvec), - ywfvec(other.ywfvec), - zwfvec(other.zwfvec), - fNPXem(other.fNPXem), - fNPYem(other.fNPYem), - fNPZem(other.fNPZem), - xemvec(other.xemvec), - yemvec(other.yemvec), - zemvec(other.zemvec), - fNPXhm(other.fNPXhm), - fNPYhm(other.fNPYhm), - fNPZhm(other.fNPZhm), - xhmvec(other.xhmvec), - yhmvec(other.yhmvec), - zhmvec(other.zhmvec), - xvecef(other.xvecef), - yvecef(other.yvecef), - zvecef(other.zvecef), - vecem(other.vecem), - vechm(other.vechm), - xvecwf(other.xvecwf), - yvecwf(other.yvecwf), - zvecwf(other.zvecwf) - { } - - __hydra_dual__ - ApplyRamo_multi operator=( ApplyRamo_multi const& other){ - if(this == &other) return *this; - - fN=other.fN; - fStep=other.fStep; - fNPXef=other.fNPXef; - fNPYef=other.fNPYef; - fNPZef=other.fNPZef; - xefvec=other.xefvec; - yefvec=other.yefvec; - zefvec=other.zefvec; - fNPXwf=other.fNPXwf; - fNPYwf=other.fNPYwf; - fNPZwf=other.fNPZwf; - xwfvec=other.xwfvec; - ywfvec=other.ywfvec; - zwfvec=other.zwfvec; - fNPXem=other.fNPXem; - fNPYem=other.fNPYem; - fNPZem=other.fNPZem; - xemvec=other.xemvec; - yemvec=other.yemvec; - zemvec=other.zemvec; - fNPXef=other.fNPXef; - fNPYef=other.fNPYef; - fNPZef=other.fNPZef; - xefvec=other.xefvec; - yefvec=other.yefvec; - zefvec=other.zefvec; - xvecef=other.xvecef; - yvecef=other.yvecef; - zvecef=other.zvecef; - vecem=other.vecem; - vechm=other.vechm; - xvecwf=other.xvecwf; - yvecwf=other.yvecwf; - zvecwf=other.zvecwf; - - return *this; - } - - //function call operator needs to be callable in host and device sides - // and be constant - template - __hydra_dual__ - void operator()(Particle p){ - - //check if particle went outside the volume - if(hydra::get<_tc_isin>(p)<0){ - hydra::get<_tc_curr>(p)=0.; - return; - } - - //get particle charge and initial position - double charge=hydra::get<_tc_charge>(p); - double mobility=0; - double x=hydra::get<_tc_x>(p); - double y=hydra::get<_tc_y>(p); - double z=hydra::get<_tc_z>(p); - - double vvec[3]; - size_t ind8[8]; - - if(charge<0.){ - getallindices(vvec, ind8, xemvec, fNPXem, yemvec, fNPYem, zemvec, fNPZem, x, y, z); - mobility=interpol(vvec, vecem, ind8); - } - else{ - getallindices(vvec, ind8, xhmvec, fNPXhm, yhmvec, fNPYhm, zhmvec, fNPZhm, x, y, z); - mobility=interpol(vvec, vechm, ind8); - } - - //get electric field - getallindices(vvec, ind8, xefvec, fNPXef, yefvec, fNPYef, zefvec, fNPZef, x, y, z); - double ex = interpol(vvec, xvecef, ind8); - double ey = interpol(vvec, yvecef, ind8); - double ez = interpol(vvec, zvecef, ind8); - - //get weighting field - getallindices(vvec, ind8, xwfvec, fNPXwf, ywfvec, fNPYwf, zwfvec, fNPZwf, x, y, z); - double wfx = interpol(vvec, xvecwf, ind8); - double wfy = interpol(vvec, yvecwf, ind8); - double wfz = interpol(vvec, zvecwf, ind8); - - //first step only - double k1x=charge*mobility*ex; - double k1y=charge*mobility*ey; - double k1z=charge*mobility*ez; - - //set diff velocity vector - double veldiff[3]={0,0,0}; - double velpoint[3]={k1x,k1y,k1z}; - getdiffvect(velpoint, veldiff, hydra::get<_tc_angle_1>(p), hydra::get<_tc_gauss_x>(p)); - - hydra::get<_tc_gauss_x>(p) = veldiff[0]; - hydra::get<_tc_gauss_y>(p) = veldiff[1]; - hydra::get<_tc_gauss_z>(p) = veldiff[2]; - - double Vdriftx=k1x+hydra::get<_tc_gauss_x>(p); - double Vdrifty=k1y+hydra::get<_tc_gauss_y>(p); - double Vdriftz=k1z+hydra::get<_tc_gauss_z>(p); - -// Vdrift=sqrt((Vdriftx*Vdriftx+Vdrifty*Vdrifty+Vdriftz*Vdriftz)); - -// MaxV=MAXVDRIFTE; -// if(charge>0) MaxV=MAXVDRIFTH; -// -// if(Vdrift>MaxV){ -// Vdriftx*=(MaxV/(Vdrift)); -// Vdrifty*=(MaxV/(Vdrift)); -// Vdriftz*=(MaxV/(Vdrift)); -// } - - //calculate current at this instant - hydra::get<5>(p)=charge*HCHARGE*((Vdriftx)*wfx+(Vdrifty)*wfy+(Vdriftz)*wfz); - - } - - - size_t fN; - double fStep; - size_t fNPXef = 0; - size_t fNPYef = 0; - size_t fNPZef = 0; - size_t fNPXwf = 0; - size_t fNPYwf = 0; - size_t fNPZwf = 0; - size_t fNPXem = 0; - size_t fNPYem = 0; - size_t fNPZem = 0; - size_t fNPXhm = 0; - size_t fNPYhm = 0; - size_t fNPZhm = 0; - double *xvecwf = nullptr; - double *yvecwf = nullptr; - double *zvecwf = nullptr; - double *xvecef = nullptr; - double *yvecef = nullptr; - double *zvecef = nullptr; - double *vecem = nullptr; - double *vechm = nullptr; - double *xefvec = nullptr; - double *yefvec = nullptr; - double *zefvec = nullptr; - double *xwfvec = nullptr; - double *ywfvec = nullptr; - double *zwfvec = nullptr; - double *xemvec = nullptr; - double *yemvec = nullptr; - double *zemvec = nullptr; - double *xhmvec = nullptr; - double *yhmvec = nullptr; - double *zhmvec = nullptr; - }; - - struct Evolve - { - - Evolve()=delete; // avoid creating objects with undefined data - - // No __hydra__dual__ here. Only allow objects to be constructed in host side. - Evolve(size_t n, double timestep, double *in_x, size_t dim_x,double *in_y, size_t dim_y,double *in_z, size_t dim_z, double *in_xef,double *in_yef, double *in_zef, double *in_em): - fN(n), - fStep(timestep), - fNPX(dim_x), - fNPY(dim_y), - fNPZ(dim_z), - xvec(in_x), - yvec(in_y), - zvec(in_z), - xvecef(in_xef), - yvecef(in_yef), - zvecef(in_zef), - vecem(in_em), - vechm(in_hm) - {} - - // No __hydra__dual__ here. Only allow objects to be constructed in host side. - Evolve(size_t n, double timestep, VecDev_t &in_x, VecDev_t &in_y, VecDev_t &in_z, VecDev_t &in_xef, VecDev_t &in_yef, VecDev_t &in_zef, VecDev_t &in_em, VecDev_t &in_hm): - fN(n), - fStep(timestep), - fNPX(in_x.size()), - fNPY(in_y.size()), - fNPZ(in_z.size()), - xvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_x.data())), - yvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_y.data())), - zvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_z.data())), - xvecef(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_xef.data())), - yvecef(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_yef.data())), - zvecef(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_zef.data())), - vecem(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_em.data())), - vechm(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_hm.data())) - {} - - // Copy constructor needs be __hydra_dual__ to be copyable - // on device and host sides, in the time of thread distribution. - // hydra, as well as thrust and std pass arguments by value. - // Obs.: RooFit passes by reference... (crazy design choice!) - __hydra_dual__ - Evolve( Evolve const& other): - fN(other.fN), - fStep(other.fStep), - fNPX(other.fNPX), - fNPY(other.fNPY), - fNPZ(other.fNPZ), - xvec(other.xvec), - yvec(other.yvec), - zvec(other.zvec), - xvecef(other.xvecef), - yvecef(other.yvecef), - zvecef(other.zvecef), - vecem(other.vecem), - vechm(other.vechm) - { } - - - __hydra_dual__ - Evolve operator=( Evolve const& other){ - if(this == &other) return *this; - - fN=other.fN; - fStep=other.fStep; - fNPX=other.fNPX; - fNPY=other.fNPY; - fNPZ=other.fNPZ; - xvec=other.xvec; - yvec=other.yvec; - zvec=other.zvec; - xvecef=other.xvecef; - yvecef=other.yvecef; - zvecef=other.zvecef; - vecem=other.vecem; - vechm=other.vechm; - return *this; - } - - //function call operator needs to be callable in host and device sides - // and be constant - template - __hydra_dual__ - void operator()(Particle p){ - - //check if particle went outside the volume - if(hydra::get<_tc_isin>(p)<0) return; - - //get particle charge and initial position - double charge=hydra::get<_tc_charge>(p); - double mobility=0; - double x=hydra::get<_tc_x>(p); - double y=hydra::get<_tc_y>(p); - double z=hydra::get<_tc_z>(p); - - //get diffusion - double vdiffx = hydra::get<_tc_gauss_x>(p); - double vdiffy = hydra::get<_tc_gauss_y>(p); - double vdiffz = hydra::get<_tc_gauss_z>(p); - - double vvec[3]; - size_t ind8[8]; - - getallindices(vvec, ind8, xvec, fNPX, yvec, fNPY, zvec, fNPZ, x, y, z); - if(charge<0.) mobility=interpol(vvec, vecem, ind8); - else mobility=interpol(vvec, vechm, ind8); - - //get electric field - double ex = interpol(vvec, xvecef, ind8); - double ey = interpol(vvec, yvecef, ind8); - double ez = interpol(vvec, zvecef, ind8); - - //first step - double k1x=charge*mobility*ex; - double k1y=charge*mobility*ey; - double k1z=charge*mobility*ez; - - //get third mobility - x = x+k1x*fStep*0.25; - y = y+k1x*fStep*0.25; - z = z+k1x*fStep*0.25; - - //get indices and versor - getallindices(vvec, ind8, xvec, fNPX, yvec, fNPY, zvec, fNPZ, x, y, z); - if(charge<0.) mobility=interpol(vvec, vecem, ind8); - else mobility=interpol(vvec, vechm, ind8); - - //get electric field - ex = interpol(vvec, xvecef, ind8); - ey = interpol(vvec, yvecef, ind8); - ez = interpol(vvec, zvecef, ind8); - - //second step - double k2x=charge*mobility*ex; - double k2y=charge*mobility*ey; - double k2z=charge*mobility*ez; - - x = x+k2x*fStep*0.25; - y = y+k2x*fStep*0.25; - z = z+k2x*fStep*0.25; - - //get indices and versor - getallindices(vvec, ind8, xvec, fNPX, yvec, fNPY, zvec, fNPZ, x, y, z); - if(charge<0.) mobility=interpol(vvec, vecem, ind8); - else mobility=interpol(vvec, vechm, ind8); - - //get electric field - ex = interpol(vvec, xvecef, ind8); - ey = interpol(vvec, yvecef, ind8); - ez = interpol(vvec, zvecef, ind8); - - //third step - double k3x=charge*mobility*ex; - double k3y=charge*mobility*ey; - double k3z=charge*mobility*ez; - - - x = x+k3x*fStep; - y = y+k3x*fStep; - z = z+k3x*fStep; - - //get indices and versor - getallindices(vvec, ind8, xvec, fNPX, yvec, fNPY, zvec, fNPZ, x, y, z); - if(charge<0.) mobility=interpol(vvec, vecem, ind8); - else mobility=interpol(vvec, vechm, ind8); - - //get electric field - ex = interpol(vvec, xvecef, ind8); - ey = interpol(vvec, yvecef, ind8); - ez = interpol(vvec, zvecef, ind8); - - //fourth step - double k4x=charge*mobility*ex; - double k4y=charge*mobility*ey; - double k4z=charge*mobility*ez; - - - //evolution - double Vdriftx=0.1666667*(k1x+2*k2x+2*k3x+k4x)+vdiffx; - double Vdrifty=0.1666667*(k1y+2*k2y+2*k3y+k4y)+vdiffy; - double Vdriftz=0.1666667*(k1z+2*k2z+2*k3z+k4z)+vdiffz; - -// Vdrift=sqrt((Vdriftx*Vdriftx+Vdrifty*Vdrifty+Vdriftz*Vdriftz)); - -// MaxV=MAXVDRIFTE; -// if(charge>0) MaxV=MAXVDRIFTH; -// -// if(Vdrift>MaxV){ -// Vdriftx*=(MaxV/(Vdrift)); -// Vdrifty*=(MaxV/(Vdrift)); -// Vdriftz*=(MaxV/(Vdrift)); -// } - - - if(xxvec[fNPX-1] || yyvec[fNPY-1] || zzvec[fNPZ-1]){ - hydra::get<_tc_isin>(p)=-1; - } - - hydra::get<_tc_x>(p)=x+Vdriftx*fStep; - hydra::get<_tc_y>(p)=y+Vdrifty*fStep; - hydra::get<_tc_z>(p)=z+Vdriftz*fStep; - - } - - - - size_t fN; - double fStep; - size_t fNPX = 0; - size_t fNPY = 0; - size_t fNPZ = 0; - double *xvecef = nullptr; - double *yvecef = nullptr; - double *zvecef = nullptr; - double *vecem = nullptr; - double *vechm = nullptr; - double *xvec = nullptr; - double *yvec = nullptr; - double *zvec = nullptr; - double *in_em = nullptr; - double *in_hm = nullptr; - }; - - struct Evolve_multi - { - - Evolve_multi()=delete; // avoid creating objects with undefined data - - // No __hydra__dual__ here. Only allow objects to be constructed in host side. - Evolve_multi(size_t n, double timestep, VecDev_t &inef_x, VecDev_t &inef_y, VecDev_t &inef_z, VecDev_t &in_xef, VecDev_t &in_yef, VecDev_t &in_zef, VecDev_t &inem_x, VecDev_t &inem_y, VecDev_t &inem_z, VecDev_t &in_em, VecDev_t &inhm_x, VecDev_t &inhm_y, VecDev_t &inhm_z, VecDev_t &in_hm, double _xmin, double _xmax, double _ymin, double _ymax, double _zmin, double _zmax): - fN(n), - fStep(timestep), - xmin(_xmin), - xmax(_xmax), - ymin(_ymin), - ymax(_ymax), - zmin(_zmin), - zmax(_zmax), - fNPXef(inef_x.size()), - fNPYef(inef_y.size()), - fNPZef(inef_z.size()), - xefvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inef_x.data())), - yefvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inef_y.data())), - zefvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inef_z.data())), - fNPXem(inem_x.size()), - fNPYem(inem_y.size()), - fNPZem(inem_z.size()), - xemvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inem_x.data())), - yemvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inem_y.data())), - zemvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inem_z.data())), - fNPXhm(inhm_x.size()), - fNPYhm(inhm_y.size()), - fNPZhm(inhm_z.size()), - xhmvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inhm_x.data())), - yhmvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inhm_y.data())), - zhmvec(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(inhm_z.data())), - xvecef(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_xef.data())), - yvecef(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_yef.data())), - zvecef(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_zef.data())), - vecem(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_em.data())), - vechm(HYDRA_EXTERNAL_NS::thrust::raw_pointer_cast(in_hm.data())) - {} - - // Copy constructor needs be __hydra_dual__ to be copyable - // on device and host sides, in the time of thread distribution. - // hydra, as well as thrust and std pass arguments by value. - // Obs.: RooFit passes by reference... (crazy design choice!) - __hydra_dual__ - Evolve_multi( Evolve_multi const& other): - fN(other.fN), - fStep(other.fStep), - xmin(other.xmin), - xmax(other.xmax), - ymin(other.ymin), - ymax(other.ymax), - zmin(other.zmin), - zmax(other.zmax), - fNPXef(other.fNPXef), - fNPYef(other.fNPYef), - fNPZef(other.fNPZef), - xefvec(other.xefvec), - yefvec(other.yefvec), - zefvec(other.zefvec), - fNPXem(other.fNPXem), - fNPYem(other.fNPYem), - fNPZem(other.fNPZem), - xemvec(other.xemvec), - yemvec(other.yemvec), - zemvec(other.zemvec), - fNPXhm(other.fNPXhm), - fNPYhm(other.fNPYhm), - fNPZhm(other.fNPZhm), - xhmvec(other.xhmvec), - yhmvec(other.yhmvec), - zhmvec(other.zhmvec), - xvecef(other.xvecef), - yvecef(other.yvecef), - zvecef(other.zvecef), - vecem(other.vecem), - vechm(other.vechm) - { } - - - __hydra_dual__ - Evolve_multi operator=( Evolve_multi const& other){ - if(this == &other) return *this; - - fN=other.fN; - fStep=other.fStep; - xmin=other.xmin; - xmax=other.xmax; - ymin=other.ymin; - ymax=other.ymax; - zmin=other.zmin; - zmax=other.zmax; - fNPXef=other.fNPXef; - fNPYef=other.fNPYef; - fNPZef=other.fNPZef; - xefvec=other.xefvec; - yefvec=other.yefvec; - zefvec=other.zefvec; - fNPXem=other.fNPXem; - fNPYem=other.fNPYem; - fNPZem=other.fNPZem; - xemvec=other.xemvec; - yemvec=other.yemvec; - zemvec=other.zemvec; - fNPXef=other.fNPXef; - fNPYef=other.fNPYef; - fNPZef=other.fNPZef; - xefvec=other.xefvec; - yefvec=other.yefvec; - zefvec=other.zefvec; - xvecef=other.xvecef; - yvecef=other.yvecef; - zvecef=other.zvecef; - vecem=other.vecem; - vechm=other.vechm; - return *this; - } - - //function call operator needs to be callable in host and device sides - // and be constant - template - __hydra_dual__ - void operator()(Particle p){ - - //check if particle went outside the volume - if(hydra::get<_tc_isin>(p)<0) return; - - //get particle charge and initial position - double charge=hydra::get<_tc_charge>(p); - double mobility=0; - double x=hydra::get<_tc_x>(p); - double y=hydra::get<_tc_y>(p); - double z=hydra::get<_tc_z>(p); - -// std::cout << p << std::endl; - - //get diffusion - double vdiffx = hydra::get<_tc_gauss_x>(p); - double vdiffy = hydra::get<_tc_gauss_y>(p); - double vdiffz = hydra::get<_tc_gauss_z>(p); - - double vvec[3]; - size_t ind8[8]; - - if(charge<0.){ - getallindices(vvec, ind8, xemvec, fNPXem, yemvec, fNPYem, zemvec, fNPZem, x, y, z); - mobility=interpol(vvec, vecem, ind8); - } - else{ - getallindices(vvec, ind8, xhmvec, fNPXhm, yhmvec, fNPYhm, zhmvec, fNPZhm, x, y, z); - mobility=interpol(vvec, vechm, ind8); - } - - //get electric field - getallindices(vvec, ind8, xefvec, fNPXef, yefvec, fNPYef, zefvec, fNPZef, x, y, z); - double ex = interpol(vvec, xvecef, ind8); - double ey = interpol(vvec, yvecef, ind8); - double ez = interpol(vvec, zvecef, ind8); - - //first step - double k1x=charge*mobility*ex; - double k1y=charge*mobility*ey; - double k1z=charge*mobility*ez; - - //get third mobility - x = x+k1x*fStep*0.25; - y = y+k1x*fStep*0.25; - z = z+k1x*fStep*0.25; - - //get indices and versor - if(charge<0.){ - getallindices(vvec, ind8, xemvec, fNPXem, yemvec, fNPYem, zemvec, fNPZem, x, y, z); - mobility=interpol(vvec, vecem, ind8); - } - else{ - getallindices(vvec, ind8, xhmvec, fNPXhm, yhmvec, fNPYhm, zhmvec, fNPZhm, x, y, z); - mobility=interpol(vvec, vechm, ind8); - } - - //get electric field - getallindices(vvec, ind8, xefvec, fNPXef, yefvec, fNPYef, zefvec, fNPZef, x, y, z); - ex = interpol(vvec, xvecef, ind8); - ey = interpol(vvec, yvecef, ind8); - ez = interpol(vvec, zvecef, ind8); - - //second step - double k2x=charge*mobility*ex; - double k2y=charge*mobility*ey; - double k2z=charge*mobility*ez; - - x = x+k2x*fStep*0.25; - y = y+k2x*fStep*0.25; - z = z+k2x*fStep*0.25; - - //get indices and versor - if(charge<0.){ - getallindices(vvec, ind8, xemvec, fNPXem, yemvec, fNPYem, zemvec, fNPZem, x, y, z); - mobility=interpol(vvec, vecem, ind8); - } - else{ - getallindices(vvec, ind8, xhmvec, fNPXhm, yhmvec, fNPYhm, zhmvec, fNPZhm, x, y, z); - mobility=interpol(vvec, vechm, ind8); - } - - //get electric field - getallindices(vvec, ind8, xefvec, fNPXef, yefvec, fNPYef, zefvec, fNPZef, x, y, z); - ex = interpol(vvec, xvecef, ind8); - ey = interpol(vvec, yvecef, ind8); - ez = interpol(vvec, zvecef, ind8); - - //third step - double k3x=charge*mobility*ex; - double k3y=charge*mobility*ey; - double k3z=charge*mobility*ez; - - - x = x+k3x*fStep; - y = y+k3x*fStep; - z = z+k3x*fStep; - - //get indices and versor - if(charge<0.){ - getallindices(vvec, ind8, xemvec, fNPXem, yemvec, fNPYem, zemvec, fNPZem, x, y, z); - mobility=interpol(vvec, vecem, ind8); - } - else{ - getallindices(vvec, ind8, xhmvec, fNPXhm, yhmvec, fNPYhm, zhmvec, fNPZhm, x, y, z); - mobility=interpol(vvec, vechm, ind8); - } - - //get electric field - getallindices(vvec, ind8, xefvec, fNPXef, yefvec, fNPYef, zefvec, fNPZef, x, y, z); - ex = interpol(vvec, xvecef, ind8); - ey = interpol(vvec, yvecef, ind8); - ez = interpol(vvec, zvecef, ind8); - - //fourth step - double k4x=charge*mobility*ex; - double k4y=charge*mobility*ey; - double k4z=charge*mobility*ez; - - - //evolution - double Vdriftx=0.1666667*(k1x+2*k2x+2*k3x+k4x)+vdiffx; - double Vdrifty=0.1666667*(k1y+2*k2y+2*k3y+k4y)+vdiffy; - double Vdriftz=0.1666667*(k1z+2*k2z+2*k3z+k4z)+vdiffz; - -// Vdrift=sqrt((Vdriftx*Vdriftx+Vdrifty*Vdrifty+Vdriftz*Vdriftz)); - -// MaxV=MAXVDRIFTE; -// if(charge>0) MaxV=MAXVDRIFTH; -// -// if(Vdrift>MaxV){ -// Vdriftx*=(MaxV/(Vdrift)); -// Vdrifty*=(MaxV/(Vdrift)); -// Vdriftz*=(MaxV/(Vdrift)); -// } - - - if(xxmax || yymax || zzmax){ - hydra::get<_tc_isin>(p)=-1; - } -// std::cout << Vdriftx << "\t" << Vdrifty << "\t" << Vdriftz << std::endl; - hydra::get<_tc_x>(p)=x+Vdriftx*fStep; - hydra::get<_tc_y>(p)=y+Vdrifty*fStep; - hydra::get<_tc_z>(p)=z+Vdriftz*fStep; - -// std::cout << fStep << std::endl; -// std::cout << Vdriftx << "\t" << Vdrifty << "\t" << Vdriftz << std::endl; -// std::cout << p << std::endl; -// std::cout << "gen"< +#include void signal_simulation(cfg::Setting const& root, int icfg=-1){ @@ -47,21 +48,21 @@ void signal_simulation(cfg::Setting const& root, int icfg=-1){ //3d hist of efield TH3D* h3=NULL; - VecDev_t _xvec, _yvec, _zvec; + VecDev_t _xvec, _yvec, _zvec; //efield vecs - VecDev_t _xvec_ef, _yvec_ef, _zvec_ef; - VecDev_t efieldvecx, efieldvecy, efieldvecz; + VecDev_t _xvec_ef, _yvec_ef, _zvec_ef; + VecDev_t efieldvecx, efieldvecy, efieldvecz; //wfield vecs - VecDev_t _xvec_wf, _yvec_wf, _zvec_wf; - VecDev_t wfieldvecx, wfieldvecy, wfieldvecz; + VecDev_t _xvec_wf, _yvec_wf, _zvec_wf; + VecDev_t wfieldvecx, wfieldvecy, wfieldvecz; //emob vecs - VecDev_t _xvec_emob, _yvec_emob, _zvec_emob; - VecDev_t emobvec; + VecDev_t _xvec_emob, _yvec_emob, _zvec_emob; + VecDev_t emobvec; //hmob vecs - VecDev_t _xvec_hmob, _yvec_hmob, _zvec_hmob; - VecDev_t hmobvec; + VecDev_t _xvec_hmob, _yvec_hmob, _zvec_hmob; + VecDev_t hmobvec; double _xmin=0, _xmax=0, _ymin=0, _ymax=0, _zmin=0, _zmax=0; //load physics maps @@ -76,9 +77,10 @@ void signal_simulation(cfg::Setting const& root, int icfg=-1){ //get physicsmap path std::string physmap_path=root["PhysicsMaps"]["map"]; maps::physmap pm(physmap_path); + pm.prepare(); - pm.fillspacepoints(_xvec,_yvec,_zvec); - pm.fillallquantities(efieldvecx,efieldvecy,efieldvecz,emobvec,hmobvec,wfieldvecx,wfieldvecy,wfieldvecz); + pm.fillspacepoints >(_xvec,_yvec,_zvec); + pm.fillallquantities >(efieldvecx,efieldvecy,efieldvecz,emobvec,hmobvec,wfieldvecx,wfieldvecy,wfieldvecz); _xmin = *(pm.getx().begin()); _xmax = *(pm.getx().rbegin()); @@ -89,8 +91,8 @@ void signal_simulation(cfg::Setting const& root, int icfg=-1){ for(auto cf:configlist){ if(std::get(cf.get("plot"))){ - VecHost_t efieldvecx_host, efieldvecy_host, efieldvecz_host; - pm.fillvector(efieldvecx_host,efieldvecy_host,efieldvecz_host); + VecHost_t efieldvecx_host, efieldvecy_host, efieldvecz_host; + pm.fillvector >(efieldvecx_host,efieldvecy_host,efieldvecz_host); h3=analysis::getHist3DDraw( efieldvecx_host, efieldvecy_host, efieldvecz_host, pm.getx(),pm.gety(),pm.getz()); break; } @@ -102,44 +104,51 @@ void signal_simulation(cfg::Setting const& root, int icfg=-1){ INFO_LINE("Loading separate maps") std::string efieldmap_path=root["PhysicsMaps"]["efield"]; + maps::physmap efieldm(efieldmap_path); //E field is a vector - efieldm.fillspacepoints(_xvec_ef,_yvec_ef,_zvec_ef); - efieldm.fillvector(efieldvecx,efieldvecy,efieldvecz); + efieldm.prepare(); + + efieldm.fillspacepoints >(_xvec_ef,_yvec_ef,_zvec_ef); + efieldm.fillvector >(efieldvecx,efieldvecy,efieldvecz); for(auto cf:configlist){ if(std::get(cf.get("plot"))){ - VecHost_t efieldvecx_host, efieldvecy_host, efieldvecz_host; - efieldm.fillvector(efieldvecx_host,efieldvecy_host,efieldvecz_host); + VecHost_t 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; } } + std::string wfieldmap_path=root["PhysicsMaps"]["wfield"]; maps::physmap wfieldm(wfieldmap_path); // Weighting field is a vector - wfieldm.fillspacepoints(_xvec_wf,_yvec_wf,_zvec_wf); - wfieldm.fillvector(wfieldvecx,wfieldvecy,wfieldvecz); + wfieldm.prepare(); + wfieldm.fillspacepoints >(_xvec_wf,_yvec_wf,_zvec_wf); + wfieldm.fillvector >(wfieldvecx,wfieldvecy,wfieldvecz); std::string emobmap_path=root["PhysicsMaps"]["emob"]; maps::physmap emobm(emobmap_path); // Electron mobility is a scalar - emobm.fillspacepoints(_xvec_emob,_yvec_emob,_zvec_emob); - emobm.fillscalar(emobvec); + emobm.prepare(); + emobm.fillspacepoints >(_xvec_emob,_yvec_emob,_zvec_emob); + emobm.fillscalar >(emobvec); std::string hmobmap_path=root["PhysicsMaps"]["hmob"]; maps::physmap hmobm(hmobmap_path); // Hole mobility is a scalar - hmobm.fillspacepoints(_xvec_hmob,_yvec_hmob,_zvec_hmob); - hmobm.fillscalar(hmobvec); + hmobm.prepare(); + hmobm.fillspacepoints >(_xvec_hmob,_yvec_hmob,_zvec_hmob); + hmobm.fillscalar >(hmobvec); - _xmin = std::min(*(efieldm.getx().begin()),std::min(*(wfieldm.getx().begin()),std::min(*(emobm.getx().begin()),*(hmobm.getx().begin())))); - _xmax = std::max(*(efieldm.getx().rbegin()),std::max(*(wfieldm.getx().rbegin()),std::max(*(emobm.getx().rbegin()),*(hmobm.getx().rbegin())))); - _ymin = std::min(*(efieldm.gety().begin()),std::min(*(wfieldm.gety().begin()),std::min(*(emobm.gety().begin()),*(hmobm.gety().begin())))); - _ymax = std::max(*(efieldm.gety().rbegin()),std::max(*(wfieldm.gety().rbegin()),std::max(*(emobm.gety().rbegin()),*(hmobm.gety().rbegin())))); - _zmin = std::min(*(efieldm.getz().begin()),std::min(*(wfieldm.getz().begin()),std::min(*(emobm.getz().begin()),*(hmobm.getz().begin())))); - _zmax = std::max(*(efieldm.getz().rbegin()),std::max(*(wfieldm.getz().rbegin()),std::max(*(emobm.getz().rbegin()),*(hmobm.getz().rbegin())))); + _xmin = std::max(*(efieldm.getx().begin()),std::max(*(wfieldm.getx().begin()),std::max(*(emobm.getx().begin()),*(hmobm.getx().begin())))); + _xmax = std::min(*(efieldm.getx().rbegin()),std::min(*(wfieldm.getx().rbegin()),std::min(*(emobm.getx().rbegin()),*(hmobm.getx().rbegin())))); + _ymin = std::max(*(efieldm.gety().begin()),std::max(*(wfieldm.gety().begin()),std::max(*(emobm.gety().begin()),*(hmobm.gety().begin())))); + _ymax = std::min(*(efieldm.gety().rbegin()),std::min(*(wfieldm.gety().rbegin()),std::min(*(emobm.gety().rbegin()),*(hmobm.gety().rbegin())))); + _zmin = std::max(*(efieldm.getz().begin()),std::max(*(wfieldm.getz().begin()),std::max(*(emobm.getz().begin()),*(hmobm.getz().begin())))); + _zmax = std::min(*(efieldm.getz().rbegin()),std::min(*(wfieldm.getz().rbegin()),std::min(*(emobm.getz().rbegin()),*(hmobm.getz().rbegin())))); std::cout << MIDDLEROW << std::endl; - INFO_LINE("Volume boundaries: [ "<< _xmin << " < x < "<< _xmax << " ] micron, "<<"[ "<< _ymin << " < y < "<< _ymax <<" ] micron, "<<"[ "<< _xmin << " < x < "<< _xmax <<" ] micron") + INFO_LINE("Volume boundaries: [ "<< _xmin << " < x < "<< _xmax << " ] micron, "<<"[ "<< _ymin << " < y < "<< _ymax <<" ] micron, "<<"[ "<< _zmin << " < z < "<< _zmax <<" ] micron") } } @@ -166,7 +175,6 @@ void signal_simulation(cfg::Setting const& root, int icfg=-1){ { std::cout << MIDDLEROW << std::endl; //get configuration - auto startconf = std::chrono::high_resolution_clock::now(); auto bunchsize = std::get(cf.get("bunchsize")); auto nparticles = std::get(cf.get("nparticles")); auto nsteps = std::get(cf.get("nsteps")); @@ -256,26 +264,16 @@ void signal_simulation(cfg::Setting const& root, int icfg=-1){ Generator.SetSeed(std::rand()); Generator.Gauss(0.0, 1.0, data_d.begin(_tc_gauss_x), data_d.end(_tc_gauss_x)); + Generator.SetSeed(std::rand()); - Generator.Gauss(0.0, 1.0, data_d.begin(_tc_gauss_y), data_d.end(_tc_gauss_y)); - Generator.SetSeed(std::rand()); - Generator.Gauss(0.0, 1.0, data_d.begin(_tc_gauss_z), data_d.end(_tc_gauss_z)); - Generator.SetSeed(std::rand()); - Generator.Gauss(0.0, 2*PI, data_d.begin(_tc_angle_1), data_d.end(_tc_angle_1)); + Generator.Uniform(0.0, 2*PI, data_d.begin(_tc_angle_1), data_d.end(_tc_angle_1)); Generator.SetSeed(std::rand()); - Generator.Gauss(0.0, 2*PI, data_d.begin(_tc_angle_2), data_d.end(_tc_angle_2)); - - - hydra::for_each(data_d, evolve::SetVeldiff_multi(niter, - timestep, - temperature_nmob, - _xvec_emob, _yvec_emob, _zvec_emob, - _xvec_hmob, _yvec_hmob, _zvec_hmob, - emobvec,hmobvec)); + Generator.Uniform(0.0, PI, data_d.begin(_tc_angle_2), data_d.end(_tc_angle_2)); hydra::for_each(data_d, evolve::ApplyRamo_multi(niter, timestep, + temperature_nmob, _xvec_ef, _yvec_ef, _zvec_ef, efieldvecx, efieldvecy, efieldvecz, _xvec_emob, _yvec_emob, _zvec_emob, emobvec, _xvec_hmob, _yvec_hmob, _zvec_hmob, hmobvec, @@ -315,25 +313,16 @@ void signal_simulation(cfg::Setting const& root, int icfg=-1){ Generator.SetSeed(std::rand()); Generator.Gauss(0.0, 1.0, data_d.begin(_tc_gauss_x), data_d.end(_tc_gauss_x)); + Generator.SetSeed(std::rand()); - Generator.Gauss(0.0, 1.0, data_d.begin(_tc_gauss_y), data_d.end(_tc_gauss_y)); - Generator.SetSeed(std::rand()); - Generator.Gauss(0.0, 1.0, data_d.begin(_tc_gauss_z), data_d.end(_tc_gauss_z)); - Generator.SetSeed(std::rand()); - Generator.Gauss(0.0, 2*PI, data_d.begin(_tc_angle_1), data_d.end(_tc_angle_1)); + Generator.Uniform(0.0, 2*PI, data_d.begin(_tc_angle_1), data_d.end(_tc_angle_1)); Generator.SetSeed(std::rand()); - Generator.Gauss(0.0, 2*PI, data_d.begin(_tc_angle_2), data_d.end(_tc_angle_2)); - - - hydra::for_each(data_d, evolve::SetVeldiff(niter, - timestep, - temperature_nmob, - _xvec, _yvec, _zvec, - emobvec,hmobvec)); + Generator.Uniform(0.0, PI, data_d.begin(_tc_angle_2), data_d.end(_tc_angle_2)); hydra::for_each(data_d, evolve::ApplyRamo(niter, timestep, + temperature_nmob, _xvec, _yvec, _zvec, efieldvecx, efieldvecy, efieldvecz, emobvec,hmobvec, diff --git a/include/analysis.h b/include/analysis.h index 1cb350a..97a0fb1 100644 --- a/include/analysis.h +++ b/include/analysis.h @@ -1,6 +1,6 @@ /*---------------------------------------------------------------------------- * - * Copyright (C) 2018 Andrea Contu e Angelo Loi + * Copyright (C) 2018-2019 Andrea Contu e Angelo Loi * * This file is part of TCode software. * @@ -22,9 +22,12 @@ * analysis.h * * Created on: 12/11/2018 - * Author: Andrea Contu + * Author: Andrea Contu */ +#ifndef __ANALYSIS_H__ +#define __ANALYSIS_H__ + //ROOT #include "TCanvas.h" #include "TROOT.h" @@ -746,3 +749,5 @@ namespace analysis{ } } + +#endif diff --git a/include/constants.h b/include/constants.h index 9c911e1..4f054a2 100644 --- a/include/constants.h +++ b/include/constants.h @@ -1,6 +1,6 @@ /*---------------------------------------------------------------------------- * - * Copyright (C) 2018 Andrea Contu e Angelo Loi + * Copyright (C) 2018-2019 Andrea Contu e Angelo Loi * * This file is part of TCode software. * @@ -22,11 +22,22 @@ * constants.h * * Created on: 12/11/2018 - * Author: Andrea Contu + * Author: Andrea Contu */ +#ifndef __CONSTANTS_H__ +#define __CONSTANTS_H__ + +//TCoDe Units +//length: micron +//time: s +//charge: coulomb +//mass: Kg + //physics quantities #define ECHARGE -1.60217662e-19 //coulombs #define HCHARGE 1.60217662e-19 //coulombs -#define KB 1.3806e-23 //boltzmann constant in SI +#define KB 1.3806e-23 //boltzmann constant in TCoDe units #define VLIGHT 3.0e12 //speed of light + +#endif diff --git a/include/datatypes.h b/include/datatypes.h index a260f7f..b8a4367 100644 --- a/include/datatypes.h +++ b/include/datatypes.h @@ -1,6 +1,6 @@ /*---------------------------------------------------------------------------- * - * Copyright (C) 2018 Andrea Contu e Angelo Loi + * Copyright (C) 2018-2019 Andrea Contu e Angelo Loi * * This file is part of TCode software. * @@ -22,57 +22,56 @@ * datatypes.h * * Created on: 12/11/2018 - * Author: Andrea Contu + * Author: Andrea Contu */ -//defining particle states at a given time -#define _tc_charge _0 -#define _tc_x _1 -#define _tc_y _2 -#define _tc_z _3 -#define _tc_isin _4 -#define _tc_curr _5 -#define _tc_gauss_x _6 -#define _tc_gauss_y _7 -#define _tc_gauss_z _8 -#define _tc_angle_1 _9 -#define _tc_angle_2 _10 -#define _tc_issec _11 -#define RSDIM 12 -#define RSHOSTDIM 7 -#define RSINT 8 +#ifndef __DATATYPES_H__ +#define __DATATYPES_H__ -typedef hydra::multiarray RunningStateHost_t; //for testing -typedef hydra::multiarray RunningState_t; -typedef hydra::tuple RunningTuple_t; +//defining particle states columns +#define RSDIM 12 // columns in running state +#define RSHOSTDIM 7 // columns in Universe state +#define RSINT 8 // entries in currents tuple + +//particle states +using RunningStateHost_t = hydra::multiarray; //for testing +using RunningState_t = hydra::multiarray; +using RunningTuple_t = hydra::tuple; +//default initialization RunningTuple_t RunningTuple_init(-1, 0., 0., 0., 1.,0.,0.,0.,0.,0.,0.,0.); -typedef hydra::multiarray StateDev_t; //Particles state in device (e.g. GPU) -typedef hydra::multiarray StateHost_t; //Particles state in host -typedef hydra::tuple StateTuple_t; +using StateDev_t = hydra::multiarray; //Particles state in device (e.g. GPU) +using StateHost_t = hydra::multiarray; //Particles state in host +using StateTuple_t = hydra::tuple; StateTuple_t StateTuple_init(0.,0.,0.,0.,0.,0.,0.); -typedef std::vector UniverseDev_t; //All particles in all states, it has to be in the host -typedef std::vector UniverseHost_t; //All particles in all states, it has to be in the host -typedef hydra::multiarray ReducedDataDev_t; //Reduced data container -typedef hydra::tuple ReducedTuple_t; +using UniverseDev_t = std::vector; //All particles in all states, it has to be in the host +using UniverseHost_t = std::vector; //All particles in all states, it has to be in the host +using ReducedDataDev_t = hydra::multiarray; //Reduced data container +using ReducedTuple_t = hydra::tuple; ReducedTuple_t ReducedTuple_init(0.,0.,0.,0.,0.,0.,0.,0.); -typedef hydra::host::vector VecHost_t; //vector container double -typedef hydra::device::vector VecDev_t; //vector container double - +template +using VecHost_t = hydra::host::vector; //vector container double +template +using VecDev_t = hydra::device::vector; //vector container double +//currents vector +using CurrentVec_t = std::vector; -typedef std::vector CurrentVec_t; -typedef std::array FileLine_t; -typedef hydra::tuple hydra_tuple3_t; -typedef hydra::tuple hydra_tuple4_t; -typedef hydra::tuple hydra_tuple6_t; -typedef hydra::tuple hydra_tuple7_t; -typedef hydra::tuple hydra_tuple8_t; -typedef hydra::tuple hydra_tuple9_t; -typedef hydra::tuple hydra_tuple10_t; -typedef hydra::tuple hydra_tuple11_t; -typedef hydra::tuple hydra_tuple12_t; +//defining particle state columns meaning +auto _tc_charge = hydra::placeholders::_0; +auto _tc_x = hydra::placeholders::_1; +auto _tc_y = hydra::placeholders::_2; +auto _tc_z = hydra::placeholders::_3; +auto _tc_isin = hydra::placeholders::_4; +auto _tc_curr = hydra::placeholders::_5; +auto _tc_gauss_x = hydra::placeholders::_6; +auto _tc_gauss_y = hydra::placeholders::_7; +auto _tc_gauss_z = hydra::placeholders::_8; +auto _tc_angle_1 = hydra::placeholders::_9; +auto _tc_angle_2 = hydra::placeholders::_10; +auto _tc_issec = hydra::placeholders::_11; +#endif diff --git a/include/defaults.h b/include/defaults.h index 3505254..fbedbda 100644 --- a/include/defaults.h +++ b/include/defaults.h @@ -1,6 +1,6 @@ /*---------------------------------------------------------------------------- * - * Copyright (C) 2018 Andrea Contu e Angelo Loi + * Copyright (C) 2018-2019 Andrea Contu e Angelo Loi * * This file is part of TCode software. * @@ -22,9 +22,12 @@ * defaults.h * * Created on: 12/11/2018 - * Author: Andrea Contu + * Author: Andrea Contu */ +#ifndef __DEFAULTS_H__ +#define __DEFAULTS_H__ + //project name #define PROJECT_NAME "TCoDe" #define VERSION "0.1-alpha" @@ -62,3 +65,4 @@ #define DEFAULT_Y 0. //default deposit shift in y direction #define DEFAULT_Z 0. //default deposit shift in z direction +#endif diff --git a/include/define.h b/include/define.h deleted file mode 100644 index 0ba3782..0000000 --- a/include/define.h +++ /dev/null @@ -1,81 +0,0 @@ -#ifndef __MULTI_MAP__ - -#define NPXEF NPX -#define NPYEF NPY -#define NPZEF NPZ -#define XDIMEF XDIM -#define YDIMEF YDIM -#define ZDIMEF ZDIM -#define XVECEF XVEC -#define YVECEF YVEC -#define ZVECEF ZVEC -#define XMINEF XMIN -#define YMINEF YMIN -#define ZMINEF ZMIN -#define XMAXEF XMAX -#define YMAXEF YMAX -#define ZMAXEF ZMAX -#define NPOINTSEF NPOINTS -#define PHYSMAPEF PHYSMAP - -#define NPXEM NPX -#define NPYEM NPY -#define NPZEM NPZ -#define XDIMEM XDIM -#define YDIMEM YDIM -#define ZDIMEM ZDIM -#define XVECEM XVEC -#define YVECEM YVEC -#define ZVECEM ZVEC -#define XMINEM XMIN -#define YMINEM YMIN -#define ZMINEM ZMIN -#define XMAXEM XMAX -#define YMAXEM YMAX -#define ZMAXEM ZMAX -#define NPOINTSEM NPOINTS -#define PHYSMAPEM PHYSMAP - -#define NPXHM NPX -#define NPYHM NPY -#define NPZHM NPZ -#define XDIMHM XDIM -#define YDIMHM YDIM -#define ZDIMHM ZDIM -#define XVECHM XVEC -#define YVECHM YVEC -#define ZVECHM ZVEC -#define XMINHM XMIN -#define YMINHM YMIN -#define ZMINHM ZMIN -#define XMAXHM XMAX -#define YMAXHM YMAX -#define ZMAXHM ZMAX -#define NPOINTSHM NPOINTS -#define PHYSMAPHM PHYSMAP - -#define NPXWF NPX -#define NPYWF NPY -#define NPZWF NPZ -#define XDIMWF XDIM -#define YDIMWF YDIM -#define ZDIMWF ZDIM -#define XVECWF XVEC -#define YVECWF YVEC -#define ZVECWF ZVEC -#define XMINWF XMIN -#define YMINWF YMIN -#define ZMINWF ZMIN -#define XMAXWF XMAX -#define YMAXWF YMAX -#define ZMAXWF ZMAX -#define NPOINTSWF NPOINTS -#define PHYSMAPWF PHYSMAP - -#define VECSDIM 3*2 - -#else - -#define VECSDIM 12*2 - -#endif diff --git a/include/deposit.h b/include/deposit.h new file mode 100644 index 0000000..85c04df --- /dev/null +++ b/include/deposit.h @@ -0,0 +1,133 @@ +/*---------------------------------------------------------------------------- + * + * Copyright (C) 2018-2019 Andrea Contu e Angelo Loi + * + * This file is part of TCode software. + * + * TCode is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * TCode is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with TCode. If not, see . + * + *---------------------------------------------------------------------------*/ +/* + * deposit.h + * + * Created on: 13/01/2019 + * Author: Andrea Contu + */ + +#ifndef __DEPOSIT_H__ +#define __DEPOSIT_H__ + +namespace deposit{ + + //class implementing deposit on host side + + class owndeposit{ + private: + //static are useful not to reload deposits if not needed + static std::string _static_path; + static RunningStateHost_t _static_data; + static size_t _static_nparticles; + static double _static_length; + + std::string _path; + RunningStateHost_t _data; + size_t _nparticles; + double _length=1.; + + public: + owndeposit(std::string path, size_t nparticles, double length){ + this->init(path, nparticles, length); + } + + owndeposit(singleconf cf){ + this->init(std::get(cf.get("path")), std::get(cf.get("nparticles")),std::get(cf.get("length"))); + } + //initialise + void init(std::string path, size_t nparticles, double length){ + if(_path==DEFAULT_PATH){ + if(nparticles!=_static_nparticles || length!=_static_length){ + _data = this->gendummy(nparticles,length); + } + else{ + _data.resize(_nparticles); + hydra::copy(_static_data,_data); + } + _nparticles = nparticles; + _length = length; + } + else{ + if(path!=_static_path){ + _data = loaddata::loadfile(path); + _nparticles = _data.size(); + _static_data.resize(_nparticles); + hydra::copy(_data,_static_data); + } + else{ + _data.resize(_nparticles); + hydra::copy(_static_data,_data); + _nparticles = _data.size(); + } + } + } + //get length + double getlength(){return _length;} + + //get nparticles + double getnparticles(){return _nparticles;} + + //generate dummy deposit + RunningStateHost_t gendummy(size_t nparticles, double length){ + INFO_LINE("Generating 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; + } + size_t halfsize=np/2; + RunningStateHost_t data_d(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.)); + + + //distribute them randomly in a line along Z + if(length>0.){ + int seed=0; + Generator.SetSeed(seed); // IMPORTANT, otherwise the seed stays the same + Generator.Uniform(0., length, data_de.begin(_tc_z), data_de.end(_tc_z)); + 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())); + + return data_d; + } + + //get deposit + RunningStateHost_t* getdata(){return &_data;} + + //get deposit and change if different + RunningStateHost_t* getdata(std::string newfile, size_t nparticles, double length=1.){ + if(_path!=newfile){ + _path=newfile; + if(newfile==DEFAULT_PATH) _data = loaddata::getDummy(nparticles,length); + else _data=loaddata::loadfile(newfile); + } + return &_data; + } + }; +} + +#endif diff --git a/include/loadfiles.h b/include/loadfiles.h index 3d745a8..bd2e73f 100644 --- a/include/loadfiles.h +++ b/include/loadfiles.h @@ -1,6 +1,6 @@ /*---------------------------------------------------------------------------- * - * Copyright (C) 2018 Andrea Contu e Angelo Loi + * Copyright (C) 2018-2019 Andrea Contu e Angelo Loi * * This file is part of TCode software. * diff --git a/include/main.inl b/include/main.inl index 5dd3bc8..4d5c192 100644 --- a/include/main.inl +++ b/include/main.inl @@ -1,6 +1,6 @@ /*---------------------------------------------------------------------------- * - * Copyright (C) 2018 Andrea Contu e Angelo Loi + * Copyright (C) 2018-2019 Andrea Contu e Angelo Loi * * This file is part of TCode software. * diff --git a/include/main_maps.inl b/include/main_maps.inl index 9d448c7..dfec02c 100644 --- a/include/main_maps.inl +++ b/include/main_maps.inl @@ -1,6 +1,6 @@ /*---------------------------------------------------------------------------- * - * Copyright (C) 2018 Andrea Contu e Angelo Loi + * Copyright (C) 2018-2019 Andrea Contu e Angelo Loi * * This file is part of TCode software. * @@ -188,41 +188,46 @@ int main(int argv, char** argc){ root.add("OutputDirectory", cfg::Setting::TypeString) = outputdir; if(PHYSMAPfile!="none"){ - root.add("pm", cfg::Setting::TypeString) = PHYSMAPfile; - root.add("pm_out", cfg::Setting::TypeString) = PHYSMAPfile_out; - root.add("pm_sc", cfg::Setting::TypeInt) = physmap_sc; - root.add("pm_sr", cfg::Setting::TypeInt) = physmap_sr; - root.add("pm_scale", cfg::Setting::TypeFloat) = physmap_scale; - root.add("pm_scalespace", cfg::Setting::TypeFloat) = physmap_scalespace; + cfg::Setting & root_pm = root["map"]; + root_pm.add("path", cfg::Setting::TypeString) = PHYSMAPfile; + root_pm.add("out", cfg::Setting::TypeString) = PHYSMAPfile_out; + root_pm.add("sc", cfg::Setting::TypeInt) = physmap_sc; + root_pm.add("sr", cfg::Setting::TypeInt) = physmap_sr; + root_pm.add("scale", cfg::Setting::TypeFloat) = physmap_scale; + root_pm.add("scalespace", cfg::Setting::TypeFloat) = physmap_scalespace; } else{ - root.add("ef", cfg::Setting::TypeString) = EFfile; - root.add("ef_out", cfg::Setting::TypeString) = EFfile_out; - root.add("ef_sc", cfg::Setting::TypeInt) = efield_sc; - root.add("ef_sr", cfg::Setting::TypeInt) = efield_sr; - root.add("ef_scale", cfg::Setting::TypeFloat) = efield_scale; - root.add("ef_scalespace", cfg::Setting::TypeFloat) = efield_scalespace; - - root.add("wf", cfg::Setting::TypeString) = WFfile; - root.add("wf_out", cfg::Setting::TypeString) = WFfile_out; - root.add("wf_sc", cfg::Setting::TypeInt) = wfield_sc; - root.add("wf_sr", cfg::Setting::TypeInt) = wfield_sr; - root.add("wf_scale", cfg::Setting::TypeFloat) = wfield_scale; - root.add("wf_scalespace", cfg::Setting::TypeFloat) = wfield_scalespace; - - root.add("emob", cfg::Setting::TypeString) = EMOBfile; - root.add("emob_out", cfg::Setting::TypeString) = EMOBfile_out; - root.add("emob_sc", cfg::Setting::TypeInt) = emob_sc; - root.add("emob_sr", cfg::Setting::TypeInt) = emob_sr; - root.add("emob_scale", cfg::Setting::TypeFloat) = emob_scale; - root.add("emob_scalespace", cfg::Setting::TypeFloat) = emob_scalespace; - - root.add("hmob", cfg::Setting::TypeString) = HMOBfile; - root.add("hmob_out", cfg::Setting::TypeString) = HMOBfile_out; - root.add("hmob_sc", cfg::Setting::TypeInt) = hmob_sc; - root.add("hmob_sr", cfg::Setting::TypeInt) = hmob_sr; - root.add("hmob_scale", cfg::Setting::TypeFloat) = hmob_scale; - root.add("hmob_scalespace", cfg::Setting::TypeFloat) = hmob_scalespace; + cfg::Setting & root_efield = root["efield"]; + root_efield.add("path", cfg::Setting::TypeString) = EFfile; + root_efield.add("out", cfg::Setting::TypeString) = EFfile_out; + root_efield.add("sc", cfg::Setting::TypeInt) = efield_sc; + root_efield.add("sr", cfg::Setting::TypeInt) = efield_sr; + root_efield.add("scale", cfg::Setting::TypeFloat) = efield_scale; + root_efield.add("scalespace", cfg::Setting::TypeFloat) = efield_scalespace; + + cfg::Setting & root_wfield = root["wfield"]; + root_wfield.add("path", cfg::Setting::TypeString) = WFfile; + root_wfield.add("out", cfg::Setting::TypeString) = WFfile_out; + root_wfield.add("sc", cfg::Setting::TypeInt) = wfield_sc; + root_wfield.add("sr", cfg::Setting::TypeInt) = wfield_sr; + root_wfield.add("scale", cfg::Setting::TypeFloat) = wfield_scale; + root_wfield.add("scalespace", cfg::Setting::TypeFloat) = wfield_scalespace; + + cfg::Setting & root_emob = root["emob"]; + root_emob.add("path", cfg::Setting::TypeString) = EMOBfile; + root_emob.add("out", cfg::Setting::TypeString) = EMOBfile_out; + root_emob.add("sc", cfg::Setting::TypeInt) = emob_sc; + root_emob.add("sr", cfg::Setting::TypeInt) = emob_sr; + root_emob.add("scale", cfg::Setting::TypeFloat) = emob_scale; + root_emob.add("scalespace", cfg::Setting::TypeFloat) = emob_scalespace; + + cfg::Setting & root_hmob = root["hmob"]; + root_hmob.add("path", cfg::Setting::TypeString) = HMOBfile; + root_hmob.add("out", cfg::Setting::TypeString) = HMOBfile_out; + root_hmob.add("sc", cfg::Setting::TypeInt) = hmob_sc; + root_hmob.add("sr", cfg::Setting::TypeInt) = hmob_sr; + root_hmob.add("scale", cfg::Setting::TypeFloat) = hmob_scale; + root_hmob.add("scalespace", cfg::Setting::TypeFloat) = hmob_scalespace; } maps::preparemaps(root); diff --git a/include/preparemaps.h b/include/preparemaps.h index e9dc2bf..9534e8c 100644 --- a/include/preparemaps.h +++ b/include/preparemaps.h @@ -1,6 +1,6 @@ /*---------------------------------------------------------------------------- * - * Copyright (C) 2018 Andrea Contu e Angelo Loi + * Copyright (C) 2018-2019 Andrea Contu e Angelo Loi * * This file is part of TCode software. * @@ -22,9 +22,12 @@ * preparemaps.h * * Created on: 12/11/2018 - * Author: Andrea Contu + * Author: Andrea Contu */ -#include + +#ifndef __PREPAREMAPS_H__ +#define __PREPAREMAPS_H__ + #include using namespace hydra::placeholders; namespace maps{ @@ -38,30 +41,6 @@ namespace maps{ } }; - struct myhashing{ - myhashing()=delete; - - myhashing(double d1, double d2): - dim1(d1), - dim2(d2) - {} - - __hydra_dual__ - myhashing(myhashing const& other): - dim1(other.dim1), - dim2(other.dim2) - {} - - template - auto operator()(T t){ - return hydra::get<0>(t)+hydra::get<1>(t)*dim1+hydra::get<2>(t)*dim1*dim2; - } - - - double dim1; - double dim2; - }; - //class for single physics map template class physmap{ @@ -95,7 +74,15 @@ namespace maps{ public: + + physmap(){ + } + physmap(std::string s, size_t skipcolumns=0, size_t skiprows=0, double scale=1., double scalespace=1., bool sparse=false){ + this->init(s,skipcolumns,skiprows,scale,scalespace,sparse); + } + + void init(std::string s, size_t skipcolumns=0, size_t skiprows=0, double scale=1., double scalespace=1., bool sparse=false){ _filename=s; _cskip=skipcolumns; _rskip=skiprows; @@ -103,6 +90,9 @@ namespace maps{ _gridyfied=false; _scalefield=scale; _scalespace=scalespace; + } + + void prepare(){ this->checkfile(); this->loaddata(); if(_isown) this->calcgrid(); @@ -112,21 +102,25 @@ namespace maps{ this->calcgrid(); } } + void set_scalespace(double A){_scalespace=A;} + void set_scale(double A){_scalefield=A;} + double get_scalespace(double A){return _scalespace;} + double get_scale(double A){return _scalefield;} + void gridify(std::set xs, std::set ys, std::set zs); //get in a grid + void write(std::string outputpath); + hydra::multiarray getdata(){return _data;} std::set getx(){return _gridx;} std::set gety(){return _gridy;} std::set getz(){return _gridz;} template void fillspacepoints(Container &vx, Container &vy, Container &vz){ - vx.reserve(_gridx.size()); - vy.reserve(_gridy.size()); - vz.reserve(_gridz.size()); vx.resize(_gridx.size()); vy.resize(_gridy.size()); vz.resize(_gridz.size()); @@ -137,16 +131,12 @@ namespace maps{ template void fillscalar(Container &v){ - v.reserve(_data.size()); v.resize(_data.size()); hydra::copy(hydra::make_range(_data.begin(_3),_data.end(_3)),v); } template void fillvector(Container &vx, Container &vy, Container &vz){ - vx.reserve(_data.size()); - vy.reserve(_data.size()); - vz.reserve(_data.size()); vx.resize(_data.size()); vy.resize(_data.size()); vz.resize(_data.size()); @@ -157,21 +147,13 @@ namespace maps{ template void fillallquantities(Container &ex, Container &ey, Container &ez, Container &emob, Container &hmob, Container &wx, Container &wy, Container &wz){ - ex.reserve(_data.size()); - ey.reserve(_data.size()); - ez.reserve(_data.size()); ex.resize(_data.size()); ey.resize(_data.size()); ez.resize(_data.size()); - wx.reserve(_data.size()); - wy.reserve(_data.size()); - wz.reserve(_data.size()); wx.resize(_data.size()); wy.resize(_data.size()); wz.resize(_data.size()); - emob.reserve(_data.size()); emob.resize(_data.size()); - hmob.reserve(_data.size()); hmob.resize(_data.size()); hydra::copy(hydra::make_range(_data.begin(_3),_data.end(_3)),ex); hydra::copy(hydra::make_range(_data.begin(_4),_data.end(_4)),ey); @@ -411,6 +393,109 @@ namespace maps{ } + + //class to manage map configuration + class mapconfig{ + private: + bool _multimap = false; + std::tuple _borders; + + maps::physmap _pm; + maps::physmap _efieldm; + maps::physmap _wfieldm; + maps::physmap _emobm; + maps::physmap _hmobm; + + VecDev_t _xvec, _yvec, _zvec; + //efield vecs + VecDev_t _xvec_ef, _yvec_ef, _zvec_ef; + VecDev_t efieldvecx, efieldvecy, efieldvecz; + //wfield vecs + VecDev_t _xvec_wf, _yvec_wf, _zvec_wf; + VecDev_t wfieldvecx, wfieldvecy, wfieldvecz; + + //emob vecs + VecDev_t _xvec_emob, _yvec_emob, _zvec_emob; + VecDev_t emobvec; + + //hmob vecs + VecDev_t _xvec_hmob, _yvec_hmob, _zvec_hmob; + VecDev_t hmobvec; + + public: + mapconfig(){} + + void load(std::string map_path){ + INFO_LINE("Loading single map") + _pm.init(map_path); + _pm.prepare(); + _pm.fillspacepoints >(_xvec,_yvec,_zvec); + _pm.fillallquantities >(efieldvecx,efieldvecy,efieldvecz,emobvec,hmobvec,wfieldvecx,wfieldvecy,wfieldvecz); + std::get<0>(_borders) = *(_pm.getx().begin()); + std::get<1>(_borders) = *(_pm.getx().rbegin()); + std::get<2>(_borders) = *(_pm.gety().begin()); + std::get<3>(_borders) = *(_pm.gety().rbegin()); + std::get<4>(_borders) = *(_pm.getz().begin()); + std::get<5>(_borders) = *(_pm.getz().rbegin()); + _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") + } + + void load(std::string efield_path, std::string emob_path, std::string hmob_path, std::string wfield_path){ + INFO_LINE("Loading separate maps") + + + _efieldm.init(efield_path); //E field is a vector + _efieldm.prepare(); + _efieldm.fillspacepoints >(_xvec_ef,_yvec_ef,_zvec_ef); + _efieldm.fillvector >(efieldvecx,efieldvecy,efieldvecz); + +// for(auto cf:configlist){ +// if(std::get(cf.get("plot"))){ +// VecHost_t 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; +// } +// } + + + _wfieldm.init(wfield_path); // Weighting field is a vector + _wfieldm.prepare(); + _wfieldm.fillspacepoints >(_xvec_wf,_yvec_wf,_zvec_wf); + _wfieldm.fillvector >(wfieldvecx,wfieldvecy,wfieldvecz); + + + _emobm.init(emob_path); // Electron mobility is a scalar + _emobm.prepare(); + _emobm.fillspacepoints >(_xvec_emob,_yvec_emob,_zvec_emob); + _emobm.fillscalar >(emobvec); + + + _hmobm.init(hmob_path); // Hole mobility is a scalar + _hmobm.prepare(); + _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())))); + std::get<1>(_borders) = std::max(*(_efieldm.getx().rbegin()),std::max(*(_wfieldm.getx().rbegin()),std::max(*(_emobm.getx().rbegin()),*(_hmobm.getx().rbegin())))); + std::get<2>(_borders) = std::min(*(_efieldm.gety().begin()),std::min(*(_wfieldm.gety().begin()),std::min(*(_emobm.gety().begin()),*(_hmobm.gety().begin())))); + std::get<3>(_borders) = std::max(*(_efieldm.gety().rbegin()),std::max(*(_wfieldm.gety().rbegin()),std::max(*(_emobm.gety().rbegin()),*(_hmobm.gety().rbegin())))); + std::get<4>(_borders) = std::min(*(_efieldm.getz().begin()),std::min(*(_wfieldm.getz().begin()),std::min(*(_emobm.getz().begin()),*(_hmobm.getz().begin())))); + std::get<5>(_borders) = std::max(*(_efieldm.getz().rbegin()),std::max(*(_wfieldm.getz().rbegin()),std::max(*(_emobm.getz().rbegin()),*(_hmobm.getz().rbegin())))); + + _multimap = true; + + 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") + } + + //is multimap or not? + bool ismulti(){return _multimap;} + }; + //prepare maps void preparemaps(cfg::Setting const& root){ @@ -422,8 +507,9 @@ namespace maps{ if(root.exists("pm")){ INFO_LINE("Loading single map") - physmap pm(root["pm"], (int)root["pm_sc"], (int)root["pm_sr"], (double)root["pm_scale"], (double)root["pm_scalespace"]); - std::string out=root["pm_out"]; + physmap pm(root["map"]["path"], (int)root["map"]["sc"], (int)root["map"]["sr"], (double)root["map"]["scale"], (double)root["map"]["scalespace"]); + pm.prepare(); + std::string out=root["pm"]["out"]; pm.write(outputdir+(std::string)"/"+out); } @@ -431,26 +517,30 @@ namespace maps{ INFO_LINE("Loading separate maps") { - physmap ef(root["ef"], (int)root["ef_sc"], (int)root["ef_sr"], (double)root["ef_scale"], (double)root["ef_scalespace"]); - std::string out_ef=root["ef_out"]; + physmap ef(root["efield"]["path"], (int)root["efield"]["sc"], (int)root["efield"]["sr"], (double)root["efield"]["scale"], (double)root["efield"]["scalespace"]); + ef.prepare(); + std::string out_ef=root["efield"]["out"]; ef.write(outputdir+(std::string)"/"+out_ef); } { - physmap wf(root["wf"], (int)root["wf_sc"], (int)root["wf_sr"], (double)root["wf_scale"], (double)root["wf_scalespace"]); - std::string out_wf=root["wf_out"]; + physmap wf(root["wfield"]["path"], (int)root["wfield"]["sc"], (int)root["wfield"]["sr"], (double)root["wfield"]["scale"], (double)root["wfield"]["scalespace"]); + wf.prepare(); + std::string out_wf=root["wfield"]["out"]; wf.write(outputdir+(std::string)"/"+out_wf); } { - physmap emob(root["emob"], (int)root["emob_sc"], (int)root["emob_sr"], (double)root["emob_scale"], (double)root["emob_scalespace"]); - std::string out_emob=root["emob_out"]; + physmap emob(root["emob"]["path"], (int)root["emob"]["sc"], (int)root["emob"]["sr"], (double)root["emob"]["scale"], (double)root["emob"]["scalespace"]); + emob.prepare(); + std::string out_emob=root["emob"]["out"]; emob.write(outputdir+(std::string)"/"+out_emob); } { - physmap hmob(root["hmob"], (int)root["hmob_sc"], (int)root["hmob_sr"], (double)root["hmob_scale"], (double)root["hmob_scalespace"]); - std::string out_hmob=root["hmob_out"]; + physmap hmob(root["hmob"]["path"], (int)root["hmob"]["sc"], (int)root["hmob"]["sr"], (double)root["hmob"]["scale"], (double)root["hmob"]["scalespace"]); + hmob.prepare(); + std::string out_hmob=root["hmob"]["out"]; hmob.write(outputdir+(std::string)"/"+out_hmob); } } @@ -462,3 +552,5 @@ namespace maps{ std::cout << TOPROW << std::endl; } } + +#endif diff --git a/include/simulation.h b/include/simulation.h new file mode 100644 index 0000000..983b3f3 --- /dev/null +++ b/include/simulation.h @@ -0,0 +1,97 @@ +/*---------------------------------------------------------------------------- + * + * Copyright (C) 2018-2019 Andrea Contu e Angelo Loi + * + * This file is part of TCode software. + * + * TCode is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * TCode is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with TCode. If not, see . + * + *---------------------------------------------------------------------------*/ +/* + * simulation.h + * + * Created on: 13/01/2019 + * Author: Andrea Contu + */ + +#ifndef __SIMULATION_H__ +#define __SIMULATION_H__ + +#include +#include + +using namespace hydra::placeholders; + +namespace sim{ + //class implementing single simulation + class simulation{ + private: + size_t _bunchsize; + size_t _nparticles; + size_t _nsteps; + std::string _path; + std::string _nickname; + bool _plot; + bool _extrainfo; + double _temperature; + double _timestep; + double _xshift; + double _yshift; + double _zshift; + double _sig; + double _amp; + double _length; +// RunningStateHost_t _data_dinit; + bool _prepared; + RunningState_t _data; + ReducedDataDev_t _currents; + CurrentVec_t _tp_currs; + + public: + //constructors + simulation(singleconf cf){ + _bunchsize = std::get(cf.get("bunchsize")); + _nparticles = std::get(cf.get("nparticles")); + _nsteps = std::get(cf.get("nsteps")); + _path = std::get(cf.get("path")); + _nickname = std::get(cf.get("nickname")); + _plot = std::get(cf.get("plot")); + _extrainfo = std::get(cf.get("extrainfo")); + _temperature = std::get(cf.get("temperature")); + _timestep = std::get(cf.get("timestep")); + _xshift = std::get(cf.get("xshift")); + _yshift = std::get(cf.get("yshift")); + _zshift = std::get(cf.get("zshift")); + _sig = std::get(cf.get("sig")); + _amp = std::get(cf.get("amp")); + _length = std::get(cf.get("length")); + _prepared = false; + } + + + //prepare simulation + void prepare(){ + deposit::owndeposit mydep(_path,_nparticles,_length); + if(_nparticles==0 || _nparticles> mydep.getnparticles()) _nparticles = mydep.getnparticles(); + _data.resize(_nparticles); + hydra::copy(hydra::make_range(mydep.getdata()->begin(),mydep.getdata()->begin() + _nparticles),_data); // copy deposit to device + _prepared = true; +// hydra::for_each(_data, loaddata::Translate(xshift,yshift,zshift,_xmin,_xmax,_ymin,_ymax,_zmin,_zmax)); + } + + }; + +} + +#endif diff --git a/include/utils.h b/include/utils.h index 3e398f4..925b864 100644 --- a/include/utils.h +++ b/include/utils.h @@ -1,6 +1,6 @@ /*---------------------------------------------------------------------------- * - * Copyright (C) 2018 Andrea Contu e Angelo Loi + * Copyright (C) 2018-2019 Andrea Contu e Angelo Loi * * This file is part of TCode software. * @@ -22,7 +22,7 @@ * utils.h * * Created on: 12/11/2018 - * Author: Andrea Contu + * Author: Andrea Contu */ #ifndef UTILS_H_