From a587ad6b3ee5174f937c5658f1b93fbfc868282f Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Sun, 29 Mar 2020 14:45:43 +0200 Subject: [PATCH] some refactoring (add '_t' to vec3 and mat3) --- include/bounding_box.hh | 4 +- include/grid_fft.hh | 32 +++++------ include/mat3.hh | 59 ++++++------------- include/particle_generator.hh | 6 +- include/particle_plt.hh | 104 +++++++++++++++++----------------- include/vec3.hh | 56 +++++++++--------- src/grid_fft.cc | 2 +- src/ic_generator.cc | 6 +- src/plugins/transfer_CLASS.cc | 3 +- 9 files changed, 124 insertions(+), 148 deletions(-) diff --git a/include/bounding_box.hh b/include/bounding_box.hh index db0f481..3048c79 100644 --- a/include/bounding_box.hh +++ b/include/bounding_box.hh @@ -5,12 +5,12 @@ template struct bounding_box { - vec3 x1_, x2_; + vec3_t x1_, x2_; bounding_box(void) { } - bounding_box( const vec3& x1, const vec3& x2) + bounding_box( const vec3_t& x1, const vec3_t& x2) : x1_(x1), x2_(x2) { } diff --git a/include/grid_fft.hh b/include/grid_fft.hh index 8acc2bd..2cf5557 100644 --- a/include/grid_fft.hh +++ b/include/grid_fft.hh @@ -165,9 +165,9 @@ public: } template - vec3 get_r(const size_t i, const size_t j, const size_t k) const + vec3_t get_r(const size_t i, const size_t j, const size_t k) const { - vec3 rr; + vec3_t rr; rr[0] = real_t(i + local_0_start_) * dx_[0]; rr[1] = real_t(j) * dx_[1]; @@ -177,9 +177,9 @@ public: } template - vec3 get_unit_r(const size_t i, const size_t j, const size_t k) const + vec3_t get_unit_r(const size_t i, const size_t j, const size_t k) const { - vec3 rr; + vec3_t rr; rr[0] = real_t(i + local_0_start_) / real_t(n_[0]); rr[1] = real_t(j) / real_t(n_[1]); @@ -189,9 +189,9 @@ public: } template - vec3 get_unit_r_shifted(const size_t i, const size_t j, const size_t k, const vec3 s) const + vec3_t get_unit_r_shifted(const size_t i, const size_t j, const size_t k, const vec3_t s) const { - vec3 rr; + vec3_t rr; rr[0] = (real_t(i + local_0_start_) + s.x) / real_t(n_[0]); rr[1] = (real_t(j) + s.y) / real_t(n_[1]); @@ -200,9 +200,9 @@ public: return rr; } - vec3 get_cell_idx_3d(const size_t i, const size_t j, const size_t k) const + vec3_t get_cell_idx_3d(const size_t i, const size_t j, const size_t k) const { - return vec3({i + local_0_start_, j, k}); + return vec3_t({i + local_0_start_, j, k}); } size_t get_cell_idx_1d(const size_t i, const size_t j, const size_t k) const @@ -226,9 +226,9 @@ public: } template - vec3 get_k(const size_t i, const size_t j, const size_t k) const + vec3_t get_k(const size_t i, const size_t j, const size_t k) const { - vec3 kk; + vec3_t kk; if( bdistributed ){ auto ip = i + local_1_start_; kk[0] = (real_t(j) - real_t(j > nhalf_[0]) * n_[0]) * kfac_[0]; @@ -243,9 +243,9 @@ public: } template - vec3 get_k(const real_t i, const real_t j, const real_t k) const + vec3_t get_k(const real_t i, const real_t j, const real_t k) const { - vec3 kk; + vec3_t kk; if( bdistributed ){ auto ip = i + real_t(local_1_start_); kk[0] = (j - real_t(j > real_t(nhalf_[0])) * n_[0]) * kfac_[0]; @@ -264,9 +264,9 @@ public: return bdistributed? std::array({j,i+local_1_start_,k}) : std::array({i,j,k}); } - data_t get_cic( const vec3& v ) const{ + data_t get_cic( const vec3_t& v ) const{ // warning! this doesn't work with MPI - vec3 x({std::fmod(v.x/length_[0]+1.0,1.0)*n_[0], + vec3_t x({std::fmod(v.x/length_[0]+1.0,1.0)*n_[0], std::fmod(v.y/length_[1]+1.0,1.0)*n_[1], std::fmod(v.z/length_[2]+1.0,1.0)*n_[2] }); size_t ix = static_cast(x.x); @@ -290,7 +290,7 @@ public: return val; } - ccomplex_t get_cic_kspace( const vec3 x ) const{ + ccomplex_t get_cic_kspace( const vec3_t x ) const{ // warning! this doesn't work with MPI int ix = static_cast(std::floor(x.x)); int iy = static_cast(std::floor(x.y)); @@ -746,7 +746,7 @@ public: void Write_PDF(std::string ofname, int nbins = 1000, double scale = 1.0, double rhomin = 1e-3, double rhomax = 1e3); - void shift_field( const vec3& s, bool transform_back=true ) + void shift_field( const vec3_t& s, bool transform_back=true ) { FourierTransformForward(); apply_function_k_dep([&](auto x, auto k) -> ccomplex_t { diff --git a/include/mat3.hh b/include/mat3.hh index ac23069..6cf2689 100644 --- a/include/mat3.hh +++ b/include/mat3.hh @@ -4,7 +4,7 @@ #include template -class mat3{ +class mat3_t{ protected: std::array data_; gsl_matrix_view m_; @@ -37,38 +37,38 @@ protected: public: - mat3() + mat3_t() : bdid_alloc_gsl_(false) {} //! copy constructor - mat3( const mat3 &m) + mat3_t( const mat3_t &m) : data_(m.data_), bdid_alloc_gsl_(false) {} //! move constructor - mat3( mat3 &&m) + mat3_t( mat3_t &&m) : data_(std::move(m.data_)), bdid_alloc_gsl_(false) {} - //! construct mat3 from initializer list + //! construct mat3_t from initializer list template - mat3(E&&...e) + mat3_t(E&&...e) : data_{{std::forward(e)...}}, bdid_alloc_gsl_(false) {} - mat3& operator=(const mat3& m) noexcept{ + mat3_t& operator=(const mat3_t& m) noexcept{ data_ = m.data_; return *this; } - mat3& operator=(const mat3&& m) noexcept{ + mat3_t& operator=(const mat3_t&& m) noexcept{ data_ = std::move(m.data_); return *this; } //! destructor - ~mat3(){ + ~mat3_t(){ this->free_gsl(); } @@ -85,7 +85,7 @@ public: const T &operator()(size_t i, size_t j) const noexcept { return data_[3*i+j]; } //! in-place addition - mat3& operator+=( const mat3& rhs ) noexcept{ + mat3_t& operator+=( const mat3_t& rhs ) noexcept{ for (size_t i = 0; i < 9; ++i) { (*this)[i] += rhs[i]; } @@ -93,7 +93,7 @@ public: } //! in-place subtraction - mat3& operator-=( const mat3& rhs ) noexcept{ + mat3_t& operator-=( const mat3_t& rhs ) noexcept{ for (size_t i = 0; i < 9; ++i) { (*this)[i] -= rhs[i]; } @@ -104,20 +104,8 @@ public: for (size_t i = 0; i < 9; ++i) data_[i]=0; } - void eigen( vec3& evals, vec3& evec1, vec3& evec2, vec3& evec3 ) + void eigen( vec3_t& evals, vec3_t& evec1, vec3_t& evec2, vec3_t& evec3_t ) { - // for( auto x : data_ ){ - // std::cerr << x << " " ; - // } - // std::cerr << std::endl; - // resort into symmetrix matrix - // data_[8] = data_[5]; - // data_[7] = data_[4]; - // data_[6] = data_[2]; - // data_[5] = data_[4]; - // data_[4] = data_[3]; - // data_[3] = data_[1]; - this->init_gsl(); gsl_eigen_symmv (&m_.matrix, eval_, evec_, wsp_); @@ -127,17 +115,15 @@ public: evals[i] = gsl_vector_get( eval_, i ); evec1[i] = gsl_matrix_get( evec_, i, 0 ); evec2[i] = gsl_matrix_get( evec_, i, 1 ); - evec3[i] = gsl_matrix_get( evec_, i, 2 ); + evec3_t[i] = gsl_matrix_get( evec_, i, 2 ); } - - // std::cerr << "(" << evals[0] << " " << evals[1] << " " << evals[2] << ")" << std::endl; } }; template -constexpr const mat3 operator+(const mat3 &lhs, const mat3 &rhs) noexcept +constexpr const mat3_t operator+(const mat3_t &lhs, const mat3_t &rhs) noexcept { - mat3 result; + mat3_t result; for (size_t i = 0; i < 9; ++i) { result[i] = lhs[i] + rhs[i]; } @@ -146,9 +132,9 @@ constexpr const mat3 operator+(const mat3 &lhs, const mat3 &rhs) noexce // matrix - vector multiplication template -vec3 operator*( const mat3 &A, const vec3 &v ) noexcept +inline vec3_t operator*( const mat3_t &A, const vec3_t &v ) noexcept { - vec3 result; + vec3_t result; for( int mu=0; mu<3; ++mu ){ result[mu] = 0.0; for( int nu=0; nu<3; ++nu ){ @@ -158,14 +144,3 @@ vec3 operator*( const mat3 &A, const vec3 &v ) noexcept return result; } -// template -// vec3 operator*( const vec3 &v, const mat3 &A ) noexcept -// { -// vec3 result = 0.0; -// for( int mu=0; mu<3; ++mu ){ -// for( int nu=0; nu<3; ++nu ){ -// result[nu] += v[mu]*A(mu,nu); -// } -// } -// return result; -// } diff --git a/include/particle_generator.hh b/include/particle_generator.hh index 956ed28..57e8b0f 100644 --- a/include/particle_generator.hh +++ b/include/particle_generator.hh @@ -18,7 +18,7 @@ enum lattice{ lattice_rsc = 3, // RSC: refined simple cubic }; -const std::vector< std::vector> > lattice_shifts = +const std::vector< std::vector> > lattice_shifts = { // first shift must always be zero! (otherwise set_positions and set_velocities break) /* SC : */ {{0.0,0.0,0.0}}, @@ -27,7 +27,7 @@ const std::vector< std::vector> > lattice_shifts = /* RSC: */ {{0.0,0.0,0.0},{0.0,0.0,0.5},{0.0,0.5,0.0},{0.0,0.5,0.5},{0.5,0.0,0.0},{0.5,0.0,0.5},{0.5,0.5,0.0},{0.5,0.5,0.5}}, }; -const std::vector> second_lattice_shift = +const std::vector> second_lattice_shift = { /* SC : */ {0.5, 0.5, 0.5}, // this corresponds to CsCl lattice /* BCC: */ {0.5, 0.5, 0.0}, // is there a diatomic lattice with BCC base?!? @@ -81,7 +81,7 @@ void set_positions( container& particles, const lattice lattice_type, bool is_se for( size_t j=0; j(i,j,k,lattice_shifts[lattice_type][ishift] - + (is_second_lattice? second_lattice_shift[lattice_type] : vec3{0.,0.,0.}) ); + + (is_second_lattice? second_lattice_shift[lattice_type] : vec3_t{0.,0.,0.}) ); if( b64reals ){ particles.set_pos64( ipcount++, idim, pos[idim]*lunit + field.relem(i,j,k) ); }else{ diff --git a/include/particle_plt.hh b/include/particle_plt.hh index b0d3760..a6fc1ad 100644 --- a/include/particle_plt.hh +++ b/include/particle_plt.hh @@ -33,13 +33,13 @@ private: const real_t mapratio_, XmL_; Grid_FFT D_xx_, D_xy_, D_xz_, D_yy_, D_yz_, D_zz_; Grid_FFT grad_x_, grad_y_, grad_z_; - std::vector> vectk_; - std::vector> ico_, vecitk_; + std::vector> vectk_; + std::vector> ico_, vecitk_; bool is_even( int i ){ return (i%2)==0; } - bool is_in( int i, int j, int k, const mat3& M ){ - vec3 v({i,j,k}); + bool is_in( int i, int j, int k, const mat3_t& M ){ + vec3_t v({i,j,k}); auto vv = M * v; return is_even(vv.x)&&is_even(vv.y)&&is_even(vv.z); } @@ -54,22 +54,22 @@ private: //! === vectors, reciprocals and normals for the SC lattice === const int charge_fac_sc = 1; - const mat3 mat_bravais_sc{ + const mat3_t mat_bravais_sc{ 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, }; - const mat3 mat_reciprocal_sc{ + const mat3_t mat_reciprocal_sc{ twopi, 0.0, 0.0, 0.0, twopi, 0.0, 0.0, 0.0, twopi, }; - const mat3 mat_invrecip_sc{ + const mat3_t mat_invrecip_sc{ 2, 0, 0, 0, 2, 0, 0, 0, 2, }; - const std::vector> normals_sc{ + const std::vector> normals_sc{ {pi,0.,0.},{-pi,0.,0.}, {0.,pi,0.},{0.,-pi,0.}, {0.,0.,pi},{0.,0.,-pi}, @@ -78,22 +78,22 @@ private: //! === vectors, reciprocals and normals for the BCC lattice === const int charge_fac_bcc = 2; - const mat3 mat_bravais_bcc{ + const mat3_t mat_bravais_bcc{ 1.0, 0.0, 0.5, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, }; - const mat3 mat_reciprocal_bcc{ + const mat3_t mat_reciprocal_bcc{ twopi, 0.0, 0.0, 0.0, twopi, 0.0, -twopi, -twopi, fourpi, }; - const mat3 mat_invrecip_bcc{ + const mat3_t mat_invrecip_bcc{ 2, 0, 0, 0, 2, 0, 1, 1, 1, }; - const std::vector> normals_bcc{ + const std::vector> normals_bcc{ {0.,pi,pi},{0.,-pi,pi},{0.,pi,-pi},{0.,-pi,-pi}, {pi,0.,pi},{-pi,0.,pi},{pi,0.,-pi},{-pi,0.,-pi}, {pi,pi,0.},{-pi,pi,0.},{pi,-pi,0.},{-pi,-pi,0.} @@ -102,22 +102,22 @@ private: //! === vectors, reciprocals and normals for the FCC lattice === const int charge_fac_fcc = 4; - const mat3 mat_bravais_fcc{ + const mat3_t mat_bravais_fcc{ 0.0, 0.5, 0.0, 0.5, 0.0, 1.0, 0.5, 0.5, 0.0, }; - const mat3 mat_reciprocal_fcc{ + const mat3_t mat_reciprocal_fcc{ -fourpi, fourpi, twopi, 0.0, 0.0, twopi, fourpi, 0.0, -twopi, }; - const mat3 mat_invrecip_fcc{ + const mat3_t mat_invrecip_fcc{ 0, 1, 1, 1, 0, 1, 0, 2, 0, }; - const std::vector> normals_fcc{ + const std::vector> normals_fcc{ {twopi,0.,0.},{-twopi,0.,0.}, {0.,twopi,0.},{0.,-twopi,0.}, {0.,0.,twopi},{0.,0.,-twopi}, @@ -152,7 +152,7 @@ private: auto kronecker = []( int i, int j ) -> real_t { return (i==j)? 1.0 : 0.0; }; //! Ewald summation: short-range Green's function - auto add_greensftide_sr = [&]( mat3& D, const vec3& d ) -> void { + auto add_greensftide_sr = [&]( mat3_t& D, const vec3_t& d ) -> void { auto r = d.norm(); if( r< 1e-14 ) return; // return zero for r=0 @@ -170,7 +170,7 @@ private: }; //! Ewald summation: long-range Green's function - auto add_greensftide_lr = [&]( mat3& D, const vec3& k, const vec3& r ) -> void { + auto add_greensftide_lr = [&]( mat3_t& D, const vec3_t& k, const vec3_t& r ) -> void { real_t kmod2 = k.norm_squared(); real_t term = std::exp(-kmod2/(4*alpha2))*std::cos(k.dot(r)) / kmod2 * fft_norm; for( int mu=0; mu<3; ++mu ){ @@ -195,22 +195,22 @@ private: constexpr ptrdiff_t lnumber = 3, knumber = 3; const int numb = 1; //!< search radius when shifting vectors into FBZ - vectk_.assign(D_xx_.memsize(),vec3()); - ico_.assign(D_xx_.memsize(),vec3()); - vecitk_.assign(D_xx_.memsize(),vec3()); + vectk_.assign(D_xx_.memsize(),vec3_t()); + ico_.assign(D_xx_.memsize(),vec3_t()); + vecitk_.assign(D_xx_.memsize(),vec3_t()); #pragma omp parallel { //... temporary to hold values of the dynamical matrix - mat3 matD(0.0); + mat3_t matD(0.0); #pragma omp for for( ptrdiff_t i=0; i x_ijk({dx*real_t(i),dx*real_t(j),dx*real_t(k)}); - const vec3 ar = (mat_bravais * x_ijk).wrap_abs(); + const vec3_t x_ijk({dx*real_t(i),dx*real_t(j),dx*real_t(k)}); + const vec3_t ar = (mat_bravais * x_ijk).wrap_abs(); //... zero temporary matrix matD.zero(); @@ -219,8 +219,8 @@ private: for( ptrdiff_t ix=-lnumber; ix<=lnumber; ix++ ){ for( ptrdiff_t iy=-lnumber; iy<=lnumber; iy++ ){ for( ptrdiff_t iz=-lnumber; iz<=lnumber; iz++ ){ - const vec3 n_ijk({real_t(ix),real_t(iy),real_t(iz)}); - const vec3 dr(ar - mat_bravais * n_ijk); + const vec3_t n_ijk({real_t(ix),real_t(iy),real_t(iz)}); + const vec3_t dr(ar - mat_bravais * n_ijk); add_greensftide_sr(matD, dr); } } @@ -231,8 +231,8 @@ private: for( ptrdiff_t iy=-knumber; iy<=knumber; iy++ ){ for( ptrdiff_t iz=-knumber; iz<=knumber; iz++ ){ if(std::abs(ix)+std::abs(iy)+std::abs(iz) != 0){ - const vec3 k_ijk({real_t(ix)/nlattice,real_t(iy)/nlattice,real_t(iz)/nlattice}); - const vec3 ak( mat_reciprocal * k_ijk); + const vec3_t k_ijk({real_t(ix)/nlattice,real_t(iy)/nlattice,real_t(iz)/nlattice}); + const vec3_t ak( mat_reciprocal * k_ijk); add_greensftide_lr(matD, ak, ar ); } @@ -278,7 +278,7 @@ private: std::ofstream ofs2("test_brillouin.txt"); #endif - using map_t = std::map,size_t>; + using map_t = std::map,size_t>; map_t iimap; //!=== Make temporary copies before resorting to std. Fourier grid ========!// @@ -312,8 +312,8 @@ private: #pragma omp parallel { // thread private matrix representation - mat3 D; - vec3 eval, evec1, evec2, evec3; + mat3_t D; + vec3_t eval, evec1, evec2, evec3_t; #pragma omp for for( size_t i=0; i kv = D_xx_.get_k(i,j,k); + vec3_t kv = D_xx_.get_k(i,j,k); // put matrix elements into actual matrix D(0,0) = std::real(temp1.kelem(i,j,k)) / fft_norm12; @@ -333,12 +333,12 @@ private: D(2,2) = std::imag(temp3.kelem(i,j,k)) / fft_norm12; // compute eigenstructure of matrix - D.eigen(eval, evec1, evec2, evec3); - evec3 /= (twopi*ngrid_); + D.eigen(eval, evec1, evec2, evec3_t); + evec3_t /= (twopi*ngrid_); // now determine to which modes on the regular lattice this contributes - vec3 ar = kv / (twopi*ngrid_); - vec3 a(mat_reciprocal * ar); + vec3_t ar = kv / (twopi*ngrid_); + vec3_t a(mat_reciprocal * ar); // translate the k-vectors into the "candidate" FBZ for( int l1=-numb; l1<=numb; ++l1 ){ @@ -347,9 +347,9 @@ private: // need both halfs of Fourier space since we use real transforms for( int isign=0; isign<=1; ++isign ){ const real_t sign = 2.0*real_t(isign)-1.0; - const vec3 vshift({real_t(l1),real_t(l2),real_t(l3)}); + const vec3_t vshift({real_t(l1),real_t(l2),real_t(l3)}); - vec3 vectk = sign * a + mat_reciprocal * vshift; + vec3_t vectk = sign * a + mat_reciprocal * vshift; if( check_FBZ( normals, vectk ) ) { @@ -358,11 +358,11 @@ private: int iz = std::round(vectk.z*(ngrid_)/twopi); #pragma omp critical - {iimap.insert( std::pair,size_t>({ix,iy,iz}, D_xx_.get_idx(i,j,k)) );} + {iimap.insert( std::pair,size_t>({ix,iy,iz}, D_xx_.get_idx(i,j,k)) );} temp1.kelem(i,j,k) = ccomplex_t(eval[2],eval[1]); - temp2.kelem(i,j,k) = ccomplex_t(eval[0],evec3.x); - temp3.kelem(i,j,k) = ccomplex_t(evec3.y,evec3.z); + temp2.kelem(i,j,k) = ccomplex_t(eval[0],evec3_t.x); + temp3.kelem(i,j,k) = ccomplex_t(evec3_t.y,evec3_t.z); } }//sign } //l3 @@ -389,24 +389,24 @@ private: int ii = (int(i)>nlattice/2)? int(i)-nlattice : int(i); int jj = (int(j)>nlattice/2)? int(j)-nlattice : int(j); int kk = (int(k)>nlattice/2)? int(k)-nlattice : int(k); - vec3 kv({real_t(ii),real_t(jj),real_t(kk)}); + vec3_t kv({real_t(ii),real_t(jj),real_t(kk)}); - auto align_with_k = [&]( const vec3& v ) -> vec3{ + auto align_with_k = [&]( const vec3_t& v ) -> vec3_t{ return v*((v.dot(kv)<0.0)?-1.0:1.0); }; - vec3 v, l; + vec3_t v, l; map_t::iterator it; if( !is_in(i,j,k,mat_invrecip) ){ - auto average_lv = [&]( const auto& t1, const auto& t2, const auto& t3, vec3& v, vec3& l ) { + auto average_lv = [&]( const auto& t1, const auto& t2, const auto& t3, vec3_t& v, vec3_t& l ) { v = 0.0; l = 0.0; int count(0); auto add_lv = [&]( auto it ) -> void { auto q = it->second;++count; - l += vec3({std::real(t1.kelem(q)),std::imag(t1.kelem(q)),std::real(t2.kelem(q))}); - v += align_with_k(vec3({std::imag(t2.kelem(q)),std::real(t3.kelem(q)),std::imag(t3.kelem(q))})); + l += vec3_t({std::real(t1.kelem(q)),std::imag(t1.kelem(q)),std::real(t2.kelem(q))}); + v += align_with_k(vec3_t({std::imag(t2.kelem(q)),std::real(t3.kelem(q)),std::imag(t3.kelem(q))})); }; map_t::iterator it; if( (it = iimap.find({ii-1,jj,kk}))!=iimap.end() ){ add_lv(it); } @@ -423,8 +423,8 @@ private: }else{ if( (it = iimap.find({ii,jj,kk}))!=iimap.end() ){ auto q = it->second; - l = vec3({std::real(temp1.kelem(q)),std::imag(temp1.kelem(q)),std::real(temp2.kelem(q))}); - v = align_with_k(vec3({std::imag(temp2.kelem(q)),std::real(temp3.kelem(q)),std::imag(temp3.kelem(q))})); + l = vec3_t({std::real(temp1.kelem(q)),std::imag(temp1.kelem(q)),std::real(temp2.kelem(q))}); + v = align_with_k(vec3_t({std::imag(temp2.kelem(q)),std::real(temp3.kelem(q)),std::imag(temp3.kelem(q))})); } } D_xx_.kelem(i,j,k) = l[0]; @@ -443,13 +443,13 @@ private: for( size_t j=0; j kv = D_xx_.get_k(i,j,k); + vec3_t kv = D_xx_.get_k(i,j,k); double mu1 = std::real(D_xx_.kelem(i,j,k)); // double mu2 = std::real(D_xy_.kelem(i,j,k)); // double mu3 = std::real(D_xz_.kelem(i,j,k)); - vec3 evec1({std::real(D_yy_.kelem(i,j,k)),std::real(D_yz_.kelem(i,j,k)),std::real(D_zz_.kelem(i,j,k))}); + vec3_t evec1({std::real(D_yy_.kelem(i,j,k)),std::real(D_yz_.kelem(i,j,k)),std::real(D_zz_.kelem(i,j,k))}); evec1 /= evec1.norm(); // /////////////////////////////////// @@ -457,7 +457,7 @@ private: real_t kr = kv.norm(), kphi = kr>0.0? std::atan2(kv.y,kv.x) : 0.0, ktheta = kr>0.0? std::acos( kv.z / kr ): 0.0; real_t st = std::sin(ktheta), ct = std::cos(ktheta), sp = std::sin(kphi), cp = std::cos(kphi); - vec3 e_r( st*cp, st*sp, ct), e_theta( ct*cp, ct*sp, -st), e_phi( -sp, cp, 0.0 ); + vec3_t e_r( st*cp, st*sp, ct), e_theta( ct*cp, ct*sp, -st), e_phi( -sp, cp, 0.0 ); // re-normalise to that longitudinal amplitude is exact double renorm = evec1.dot( e_r ); if( renorm < 0.01 ) renorm = 1.0; diff --git a/include/vec3.hh b/include/vec3.hh index 4e72d81..3d1fe44 100644 --- a/include/vec3.hh +++ b/include/vec3.hh @@ -1,5 +1,5 @@ /*******************************************************************\ - vec3.hh - This file is part of MUSIC2 - + vec3_t.hh - This file is part of MUSIC2 - a code to generate initial conditions for cosmological simulations CHANGELOG (only majors, for details see repo): @@ -9,7 +9,7 @@ //! implements a simple class of 3-vectors of arbitrary scalar type template< typename T > -class vec3{ +class vec3_t{ private: //! holds the data std::array data_; @@ -19,27 +19,27 @@ public: T &x,&y,&z; //! empty constructor - vec3() + vec3_t() : x(data_[0]),y(data_[1]),z(data_[2]){} //! copy constructor - vec3( const vec3 &v) + vec3_t( const vec3_t &v) : data_(v.data_), x(data_[0]),y(data_[1]),z(data_[2]){} //! copy constructor for non-const reference, needed to avoid variadic template being called for non-const reference - vec3( vec3& v) + vec3_t( vec3_t& v) : data_(v.data_), x(data_[0]),y(data_[1]),z(data_[2]){} //! move constructor - vec3( vec3 &&v) + vec3_t( vec3_t &&v) : data_(std::move(v.data_)), x(data_[0]), y(data_[1]), z(data_[2]){} - //! construct vec3 from initializer list + //! construct vec3_t from initializer list template - vec3(E&&...e) + vec3_t(E&&...e) : data_{{std::forward(e)...}}, x{data_[0]}, y{data_[1]}, z{data_[2]} {} - // vec3( T a, T b, T c ) + // vec3_t( T a, T b, T c ) // : data_{{a,b,c}}, x(data_[0]), y(data_[1]), z(data_[2]){} //! bracket index access to vector components @@ -49,37 +49,37 @@ public: const T &operator[](size_t i) const noexcept { return data_[i]; } // assignment operator - vec3& operator=( const vec3& v ) noexcept { data_=v.data_; return *this; } + vec3_t& operator=( const vec3_t& v ) noexcept { data_=v.data_; return *this; } - //! implementation of summation of vec3 - vec3 operator+( const vec3& v ) const noexcept{ return vec3({x+v.x,y+v.y,z+v.z}); } + //! implementation of summation of vec3_t + vec3_t operator+( const vec3_t& v ) const noexcept{ return vec3_t({x+v.x,y+v.y,z+v.z}); } - //! implementation of difference of vec3 - vec3 operator-( const vec3& v ) const noexcept{ return vec3({x-v.x,y-v.y,z-v.z}); } + //! implementation of difference of vec3_t + vec3_t operator-( const vec3_t& v ) const noexcept{ return vec3_t({x-v.x,y-v.y,z-v.z}); } //! implementation of unary negative - vec3 operator-() const noexcept{ return vec3({-x,-y,-z}); } + vec3_t operator-() const noexcept{ return vec3_t({-x,-y,-z}); } //! implementation of scalar multiplication - vec3 operator*( T s ) const noexcept{ return vec3({x*s,y*s,z*s}); } + vec3_t operator*( T s ) const noexcept{ return vec3_t({x*s,y*s,z*s}); } //! implementation of scalar division - vec3 operator/( T s ) const noexcept{ return vec3({x/s,y/s,z/s}); } + vec3_t operator/( T s ) const noexcept{ return vec3_t({x/s,y/s,z/s}); } //! implementation of += operator - vec3& operator+=( const vec3& v ) noexcept{ x+=v.x; y+=v.y; z+=v.z; return *this; } + vec3_t& operator+=( const vec3_t& v ) noexcept{ x+=v.x; y+=v.y; z+=v.z; return *this; } //! implementation of -= operator - vec3& operator-=( const vec3& v ) noexcept{ x-=v.x; y-=v.y; z-=v.z; return *this; } + vec3_t& operator-=( const vec3_t& v ) noexcept{ x-=v.x; y-=v.y; z-=v.z; return *this; } //! multiply with scalar - vec3& operator*=( T s ) noexcept{ x*=s; y*=s; z*=s; return *this; } + vec3_t& operator*=( T s ) noexcept{ x*=s; y*=s; z*=s; return *this; } //! divide by scalar - vec3& operator/=( T s ) noexcept{ x/=s; y/=s; z/=s; return *this; } + vec3_t& operator/=( T s ) noexcept{ x/=s; y/=s; z/=s; return *this; } //! compute dot product with another vector - T dot(const vec3 &a) const noexcept + T dot(const vec3_t &a) const noexcept { return data_[0] * a.data_[0] + data_[1] * a.data_[1] + data_[2] * a.data_[2]; } @@ -91,19 +91,19 @@ public: T norm(void) const noexcept { return std::sqrt( this->norm_squared() ); } //! wrap absolute vector to box of size p - vec3& wrap_abs( T p = 1.0 ) noexcept{ + vec3_t& wrap_abs( T p = 1.0 ) noexcept{ for( auto& x : data_ ) x = std::fmod( 2*p + x, p ); return *this; } //! wrap relative vector to box of size p - vec3& wrap_rel( T p = 1.0 ) noexcept{ + vec3_t& wrap_rel( T p = 1.0 ) noexcept{ for( auto& x : data_ ) x = (x<-p/2)? x+p : (x>=p/2)? x-p : x; return *this; } - //! ordering, allows 3d sorting of vec3s - bool operator<( const vec3& o ) const noexcept{ + //! ordering, allows 3d sorting of vec3_ts + bool operator<( const vec3_t& o ) const noexcept{ if( x!=o.x ) return x -vec3 operator*( T s, const vec3& v ){ - return vec3({v.x*s,v.y*s,v.z*s}); +vec3_t operator*( T s, const vec3_t& v ){ + return vec3_t({v.x*s,v.y*s,v.z*s}); } diff --git a/src/grid_fft.cc b/src/grid_fft.cc index 5ae6b24..a1b1912 100644 --- a/src/grid_fft.cc +++ b/src/grid_fft.cc @@ -860,7 +860,7 @@ void Grid_FFT::Compute_PowerSpectrum(std::vector &b for (size_t iy = 0; iy < size(1); iy++) for (size_t iz = 0; iz < size(2); iz++) { - vec3 k3 = get_k(ix, iy, iz); + vec3_t k3 = get_k(ix, iy, iz); double k = k3.norm(); int idx2 = k / dk; //int((1.0f / dklog * std::log10(k / kmin))); auto z = this->kelem(ix, iy, iz); diff --git a/src/ic_generator.cc b/src/ic_generator.cc index ce86444..708e12b 100644 --- a/src/ic_generator.cc +++ b/src/ic_generator.cc @@ -320,7 +320,7 @@ int Run( ConfigFile& the_config ) if (bAddExternalTides) { - phi2.assign_function_of_grids_kdep([&](vec3 kvec, ccomplex_t pphi, ccomplex_t pphi2) { + phi2.assign_function_of_grids_kdep([&](vec3_t kvec, ccomplex_t pphi, ccomplex_t pphi2) { // sign in front of f_aniso is reversed since phi1 = -phi return pphi2 + f_aniso * (kvec[0] * kvec[0] * lss_aniso_lambda[0] + kvec[1] * kvec[1] * lss_aniso_lambda[1] + kvec[2] * kvec[2] * lss_aniso_lambda[2]) * pphi; }, @@ -569,7 +569,7 @@ int Run( ConfigFile& the_config ) + lg.gradient(idimp,tmp.get_k3(i,j,k)) * A3[idimpp]->kelem(idx) - lg.gradient(idimpp,tmp.get_k3(i,j,k)) * A3[idimp]->kelem(idx) ); if( bDoBaryons ){ - vec3 kvec = phi.get_k(i,j,k); + vec3_t kvec = phi.get_k(i,j,k); real_t k2 = kvec.norm_squared(), kmod = std::sqrt(k2); double ampldiff = ((this_species == cosmo_species::dm)? the_cosmo_calc->GetAmplitude(kmod, cdm0) : (this_species == cosmo_species::baryon)? the_cosmo_calc->GetAmplitude(kmod, baryon0) : @@ -618,7 +618,7 @@ int Run( ConfigFile& the_config ) + vfac3 * (lg.gradient(idimp,tmp.get_k3(i,j,k)) * A3[idimpp]->kelem(idx) - lg.gradient(idimpp,tmp.get_k3(i,j,k)) * A3[idimp]->kelem(idx)) ); if( bDoBaryons ){ - vec3 kvec = phi.get_k(i,j,k); + vec3_t kvec = phi.get_k(i,j,k); real_t k2 = kvec.norm_squared(), kmod = std::sqrt(k2); double ampldiff = ((this_species == cosmo_species::dm)? the_cosmo_calc->GetAmplitude(kmod, vcdm0) : (this_species == cosmo_species::baryon)? the_cosmo_calc->GetAmplitude(kmod, vbaryon0) : diff --git a/src/plugins/transfer_CLASS.cc b/src/plugins/transfer_CLASS.cc index e6e2c00..079b633 100644 --- a/src/plugins/transfer_CLASS.cc +++ b/src/plugins/transfer_CLASS.cc @@ -223,13 +223,14 @@ public: gsl_spline *splineT = nullptr; gsl_interp_accel *accT = nullptr; switch(type){ + // values at ztarget: case total: splineT = gsl_sp_dtot_; accT = gsl_ia_dtot_; break; case cdm: splineT = gsl_sp_dc_; accT = gsl_ia_dc_; break; case baryon: splineT = gsl_sp_db_; accT = gsl_ia_db_; break; case vtotal: splineT = gsl_sp_ttot_; accT = gsl_ia_ttot_; break; case vcdm: splineT = gsl_sp_tc_; accT = gsl_ia_tc_; break; case vbaryon: splineT = gsl_sp_tb_; accT = gsl_ia_tb_; break; - + // values at zstart: case total0: splineT = gsl_sp_dtot0_;accT = gsl_ia_dtot0_;break; case cdm0: splineT = gsl_sp_dc0_; accT = gsl_ia_dc0_; break; case baryon0: splineT = gsl_sp_db0_; accT = gsl_ia_db0_; break;