mirror/monofonIC
1
0
Fork 0
mirror of https://github.com/cosmo-sims/monofonIC.git synced 2024-09-19 17:03:45 +02:00
Code Releases Activity

some refactoring (add '_t' to vec3 and mat3)

Browse source
This commit is contained in:
Oliver Hahn 2020-03-29 14:45:43 +02:00
parent ab2db06990
commit a587ad6b3e
9 changed files with 124 additions and 148 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
include/bounding_box.hh
View file

@ -5,12 +5,12 @@
template <typename T> template <typename T>
struct bounding_box struct bounding_box
{ {
vec3<T> x1_, x2_; vec3_t<T> x1_, x2_;
bounding_box(void) bounding_box(void)
{ } { }
bounding_box( const vec3<T>& x1, const vec3<T>& x2) bounding_box( const vec3_t<T>& x1, const vec3_t<T>& x2)
: x1_(x1), x2_(x2) : x1_(x1), x2_(x2)
{ } { }

32
include/grid_fft.hh
View file

@ -165,9 +165,9 @@ public:
} }
template <typename ft> template <typename ft>
vec3<ft> get_r(const size_t i, const size_t j, const size_t k) const vec3_t<ft> get_r(const size_t i, const size_t j, const size_t k) const
{ {
vec3<ft> rr; vec3_t<ft> rr;
rr[0] = real_t(i + local_0_start_) * dx_[0]; rr[0] = real_t(i + local_0_start_) * dx_[0];
rr[1] = real_t(j) * dx_[1]; rr[1] = real_t(j) * dx_[1];
@ -177,9 +177,9 @@ public:
} }
template <typename ft> template <typename ft>
vec3<ft> get_unit_r(const size_t i, const size_t j, const size_t k) const vec3_t<ft> get_unit_r(const size_t i, const size_t j, const size_t k) const
{ {
vec3<ft> rr; vec3_t<ft> rr;
rr[0] = real_t(i + local_0_start_) / real_t(n_[0]); rr[0] = real_t(i + local_0_start_) / real_t(n_[0]);
rr[1] = real_t(j) / real_t(n_[1]); rr[1] = real_t(j) / real_t(n_[1]);
@ -189,9 +189,9 @@ public:
} }
template <typename ft> template <typename ft>
vec3<ft> get_unit_r_shifted(const size_t i, const size_t j, const size_t k, const vec3<real_t> s) const vec3_t<ft> get_unit_r_shifted(const size_t i, const size_t j, const size_t k, const vec3_t<real_t> s) const
{ {
vec3<ft> rr; vec3_t<ft> rr;
rr[0] = (real_t(i + local_0_start_) + s.x) / real_t(n_[0]); 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]); rr[1] = (real_t(j) + s.y) / real_t(n_[1]);
@ -200,9 +200,9 @@ public:
return rr; return rr;
} }
vec3<size_t> get_cell_idx_3d(const size_t i, const size_t j, const size_t k) const vec3_t<size_t> get_cell_idx_3d(const size_t i, const size_t j, const size_t k) const
{ {
return vec3<size_t>({i + local_0_start_, j, k}); return vec3_t<size_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 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 <typename ft> template <typename ft>
vec3<ft> get_k(const size_t i, const size_t j, const size_t k) const vec3_t<ft> get_k(const size_t i, const size_t j, const size_t k) const
{ {
vec3<ft> kk; vec3_t<ft> kk;
if( bdistributed ){ if( bdistributed ){
auto ip = i + local_1_start_; auto ip = i + local_1_start_;
kk[0] = (real_t(j) - real_t(j > nhalf_[0]) * n_[0]) * kfac_[0]; kk[0] = (real_t(j) - real_t(j > nhalf_[0]) * n_[0]) * kfac_[0];
@ -243,9 +243,9 @@ public:
} }
template <typename ft> template <typename ft>
vec3<ft> get_k(const real_t i, const real_t j, const real_t k) const vec3_t<ft> get_k(const real_t i, const real_t j, const real_t k) const
{ {
vec3<ft> kk; vec3_t<ft> kk;
if( bdistributed ){ if( bdistributed ){
auto ip = i + real_t(local_1_start_); auto ip = i + real_t(local_1_start_);
kk[0] = (j - real_t(j > real_t(nhalf_[0])) * n_[0]) * kfac_[0]; kk[0] = (j - real_t(j > real_t(nhalf_[0])) * n_[0]) * kfac_[0];
@ -264,9 +264,9 @@ public:
return bdistributed? std::array<size_t,3>({j,i+local_1_start_,k}) : std::array<size_t,3>({i,j,k}); return bdistributed? std::array<size_t,3>({j,i+local_1_start_,k}) : std::array<size_t,3>({i,j,k});
} }
data_t get_cic( const vec3<real_t>& v ) const{ data_t get_cic( const vec3_t<real_t>& v ) const{
// warning! this doesn't work with MPI // warning! this doesn't work with MPI
vec3<real_t> x({std::fmod(v.x/length_[0]+1.0,1.0)*n_[0], vec3_t<real_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.y/length_[1]+1.0,1.0)*n_[1],
std::fmod(v.z/length_[2]+1.0,1.0)*n_[2] }); std::fmod(v.z/length_[2]+1.0,1.0)*n_[2] });
size_t ix = static_cast<size_t>(x.x); size_t ix = static_cast<size_t>(x.x);
@ -290,7 +290,7 @@ public:
return val; return val;
} }
ccomplex_t get_cic_kspace( const vec3<real_t> x ) const{ ccomplex_t get_cic_kspace( const vec3_t<real_t> x ) const{
// warning! this doesn't work with MPI // warning! this doesn't work with MPI
int ix = static_cast<int>(std::floor(x.x)); int ix = static_cast<int>(std::floor(x.x));
int iy = static_cast<int>(std::floor(x.y)); int iy = static_cast<int>(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 Write_PDF(std::string ofname, int nbins = 1000, double scale = 1.0, double rhomin = 1e-3, double rhomax = 1e3);
void shift_field( const vec3<real_t>& s, bool transform_back=true ) void shift_field( const vec3_t<real_t>& s, bool transform_back=true )
{ {
FourierTransformForward(); FourierTransformForward();
apply_function_k_dep([&](auto x, auto k) -> ccomplex_t { apply_function_k_dep([&](auto x, auto k) -> ccomplex_t {

59
include/mat3.hh
View file

@ -4,7 +4,7 @@
#include <vec3.hh> #include <vec3.hh>
template<typename T> template<typename T>
class mat3{ class mat3_t{
protected: protected:
std::array<T,9> data_; std::array<T,9> data_;
gsl_matrix_view m_; gsl_matrix_view m_;
@ -37,38 +37,38 @@ protected:
public: public:
mat3() mat3_t()
: bdid_alloc_gsl_(false) : bdid_alloc_gsl_(false)
{} {}
//! copy constructor //! copy constructor
mat3( const mat3<T> &m) mat3_t( const mat3_t<T> &m)
: data_(m.data_), bdid_alloc_gsl_(false) : data_(m.data_), bdid_alloc_gsl_(false)
{} {}
//! move constructor //! move constructor
mat3( mat3<T> &&m) mat3_t( mat3_t<T> &&m)
: data_(std::move(m.data_)), bdid_alloc_gsl_(false) : data_(std::move(m.data_)), bdid_alloc_gsl_(false)
{} {}
//! construct mat3 from initializer list //! construct mat3_t from initializer list
template<typename ...E> template<typename ...E>
mat3(E&&...e) mat3_t(E&&...e)
: data_{{std::forward<E>(e)...}}, bdid_alloc_gsl_(false) : data_{{std::forward<E>(e)...}}, bdid_alloc_gsl_(false)
{} {}
mat3<T>& operator=(const mat3<T>& m) noexcept{ mat3_t<T>& operator=(const mat3_t<T>& m) noexcept{
data_ = m.data_; data_ = m.data_;
return *this; return *this;
} }
mat3<T>& operator=(const mat3<T>&& m) noexcept{ mat3_t<T>& operator=(const mat3_t<T>&& m) noexcept{
data_ = std::move(m.data_); data_ = std::move(m.data_);
return *this; return *this;
} }
//! destructor //! destructor
~mat3(){ ~mat3_t(){
this->free_gsl(); this->free_gsl();
} }
@ -85,7 +85,7 @@ public:
const T &operator()(size_t i, size_t j) const noexcept { return data_[3*i+j]; } const T &operator()(size_t i, size_t j) const noexcept { return data_[3*i+j]; }
//! in-place addition //! in-place addition
mat3<T>& operator+=( const mat3<T>& rhs ) noexcept{ mat3_t<T>& operator+=( const mat3_t<T>& rhs ) noexcept{
for (size_t i = 0; i < 9; ++i) { for (size_t i = 0; i < 9; ++i) {
(*this)[i] += rhs[i]; (*this)[i] += rhs[i];
} }
@ -93,7 +93,7 @@ public:
} }
//! in-place subtraction //! in-place subtraction
mat3<T>& operator-=( const mat3<T>& rhs ) noexcept{ mat3_t<T>& operator-=( const mat3_t<T>& rhs ) noexcept{
for (size_t i = 0; i < 9; ++i) { for (size_t i = 0; i < 9; ++i) {
(*this)[i] -= rhs[i]; (*this)[i] -= rhs[i];
} }
@ -104,20 +104,8 @@ public:
for (size_t i = 0; i < 9; ++i) data_[i]=0; for (size_t i = 0; i < 9; ++i) data_[i]=0;
} }
void eigen( vec3<T>& evals, vec3<T>& evec1, vec3<T>& evec2, vec3<T>& evec3 ) void eigen( vec3_t<T>& evals, vec3_t<T>& evec1, vec3_t<T>& evec2, vec3_t<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(); this->init_gsl();
gsl_eigen_symmv (&m_.matrix, eval_, evec_, wsp_); gsl_eigen_symmv (&m_.matrix, eval_, evec_, wsp_);
@ -127,17 +115,15 @@ public:
evals[i] = gsl_vector_get( eval_, i ); evals[i] = gsl_vector_get( eval_, i );
evec1[i] = gsl_matrix_get( evec_, i, 0 ); evec1[i] = gsl_matrix_get( evec_, i, 0 );
evec2[i] = gsl_matrix_get( evec_, i, 1 ); 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<typename T> template<typename T>
constexpr const mat3<T> operator+(const mat3<T> &lhs, const mat3<T> &rhs) noexcept constexpr const mat3_t<T> operator+(const mat3_t<T> &lhs, const mat3_t<T> &rhs) noexcept
{ {
mat3<T> result; mat3_t<T> result;
for (size_t i = 0; i < 9; ++i) { for (size_t i = 0; i < 9; ++i) {
result[i] = lhs[i] + rhs[i]; result[i] = lhs[i] + rhs[i];
} }
@ -146,9 +132,9 @@ constexpr const mat3<T> operator+(const mat3<T> &lhs, const mat3<T> &rhs) noexce
// matrix - vector multiplication // matrix - vector multiplication
template<typename T> template<typename T>
vec3<T> operator*( const mat3<T> &A, const vec3<T> &v ) noexcept inline vec3_t<T> operator*( const mat3_t<T> &A, const vec3_t<T> &v ) noexcept
{ {
vec3<T> result; vec3_t<T> result;
for( int mu=0; mu<3; ++mu ){ for( int mu=0; mu<3; ++mu ){
result[mu] = 0.0; result[mu] = 0.0;
for( int nu=0; nu<3; ++nu ){ for( int nu=0; nu<3; ++nu ){
@ -158,14 +144,3 @@ vec3<T> operator*( const mat3<T> &A, const vec3<T> &v ) noexcept
return result; return result;
} }
// template<typename T>
// vec3<T> operator*( const vec3<T> &v, const mat3<T> &A ) noexcept
// {
// vec3<T> 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;
// }

6
include/particle_generator.hh
View file

@ -18,7 +18,7 @@ enum lattice{
lattice_rsc = 3, // RSC: refined simple cubic lattice_rsc = 3, // RSC: refined simple cubic
}; };
const std::vector< std::vector<vec3<real_t>> > lattice_shifts = const std::vector< std::vector<vec3_t<real_t>> > lattice_shifts =
{ {
// first shift must always be zero! (otherwise set_positions and set_velocities break) // first shift must always be zero! (otherwise set_positions and set_velocities break)
/* SC : */ {{0.0,0.0,0.0}}, /* SC : */ {{0.0,0.0,0.0}},
@ -27,7 +27,7 @@ const std::vector< std::vector<vec3<real_t>> > 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}}, /* 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<vec3<real_t>> second_lattice_shift = const std::vector<vec3_t<real_t>> second_lattice_shift =
{ {
/* SC : */ {0.5, 0.5, 0.5}, // this corresponds to CsCl lattice /* 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?!? /* 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<field.size(1); ++j){ for( size_t j=0; j<field.size(1); ++j){
for( size_t k=0; k<field.size(2); ++k){ for( size_t k=0; k<field.size(2); ++k){
auto pos = field.template get_unit_r_shifted<real_t>(i,j,k,lattice_shifts[lattice_type][ishift] auto pos = field.template get_unit_r_shifted<real_t>(i,j,k,lattice_shifts[lattice_type][ishift]
+ (is_second_lattice? second_lattice_shift[lattice_type] : vec3<real_t>{0.,0.,0.}) ); + (is_second_lattice? second_lattice_shift[lattice_type] : vec3_t<real_t>{0.,0.,0.}) );
if( b64reals ){ if( b64reals ){
particles.set_pos64( ipcount++, idim, pos[idim]*lunit + field.relem(i,j,k) ); particles.set_pos64( ipcount++, idim, pos[idim]*lunit + field.relem(i,j,k) );
}else{ }else{

104
include/particle_plt.hh
View file

@ -33,13 +33,13 @@ private:
const real_t mapratio_, XmL_; const real_t mapratio_, XmL_;
Grid_FFT<real_t,false> D_xx_, D_xy_, D_xz_, D_yy_, D_yz_, D_zz_; Grid_FFT<real_t,false> D_xx_, D_xy_, D_xz_, D_yy_, D_yz_, D_zz_;
Grid_FFT<real_t,false> grad_x_, grad_y_, grad_z_; Grid_FFT<real_t,false> grad_x_, grad_y_, grad_z_;
std::vector<vec3<real_t>> vectk_; std::vector<vec3_t<real_t>> vectk_;
std::vector<vec3<int>> ico_, vecitk_; std::vector<vec3_t<int>> ico_, vecitk_;
bool is_even( int i ){ return (i%2)==0; } bool is_even( int i ){ return (i%2)==0; }
bool is_in( int i, int j, int k, const mat3<int>& M ){ bool is_in( int i, int j, int k, const mat3_t<int>& M ){
vec3<int> v({i,j,k}); vec3_t<int> v({i,j,k});
auto vv = M * v; auto vv = M * v;
return is_even(vv.x)&&is_even(vv.y)&&is_even(vv.z); 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 === //! === vectors, reciprocals and normals for the SC lattice ===
const int charge_fac_sc = 1; const int charge_fac_sc = 1;
const mat3<real_t> mat_bravais_sc{ const mat3_t<real_t> mat_bravais_sc{
1.0, 0.0, 0.0, 1.0, 0.0, 0.0,
0.0, 1.0, 0.0, 0.0, 1.0, 0.0,
0.0, 0.0, 1.0, 0.0, 0.0, 1.0,
}; };
const mat3<real_t> mat_reciprocal_sc{ const mat3_t<real_t> mat_reciprocal_sc{
twopi, 0.0, 0.0, twopi, 0.0, 0.0,
0.0, twopi, 0.0, 0.0, twopi, 0.0,
0.0, 0.0, twopi, 0.0, 0.0, twopi,
}; };
const mat3<int> mat_invrecip_sc{ const mat3_t<int> mat_invrecip_sc{
2, 0, 0, 2, 0, 0,
0, 2, 0, 0, 2, 0,
0, 0, 2, 0, 0, 2,
}; };
const std::vector<vec3<real_t>> normals_sc{ const std::vector<vec3_t<real_t>> normals_sc{
{pi,0.,0.},{-pi,0.,0.}, {pi,0.,0.},{-pi,0.,0.},
{0.,pi,0.},{0.,-pi,0.}, {0.,pi,0.},{0.,-pi,0.},
{0.,0.,pi},{0.,0.,-pi}, {0.,0.,pi},{0.,0.,-pi},
@ -78,22 +78,22 @@ private:
//! === vectors, reciprocals and normals for the BCC lattice === //! === vectors, reciprocals and normals for the BCC lattice ===
const int charge_fac_bcc = 2; const int charge_fac_bcc = 2;
const mat3<real_t> mat_bravais_bcc{ const mat3_t<real_t> mat_bravais_bcc{
1.0, 0.0, 0.5, 1.0, 0.0, 0.5,
0.0, 1.0, 0.5, 0.0, 1.0, 0.5,
0.0, 0.0, 0.5, 0.0, 0.0, 0.5,
}; };
const mat3<real_t> mat_reciprocal_bcc{ const mat3_t<real_t> mat_reciprocal_bcc{
twopi, 0.0, 0.0, twopi, 0.0, 0.0,
0.0, twopi, 0.0, 0.0, twopi, 0.0,
-twopi, -twopi, fourpi, -twopi, -twopi, fourpi,
}; };
const mat3<int> mat_invrecip_bcc{ const mat3_t<int> mat_invrecip_bcc{
2, 0, 0, 2, 0, 0,
0, 2, 0, 0, 2, 0,
1, 1, 1, 1, 1, 1,
}; };
const std::vector<vec3<real_t>> normals_bcc{ const std::vector<vec3_t<real_t>> normals_bcc{
{0.,pi,pi},{0.,-pi,pi},{0.,pi,-pi},{0.,-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.,-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.} {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 === //! === vectors, reciprocals and normals for the FCC lattice ===
const int charge_fac_fcc = 4; const int charge_fac_fcc = 4;
const mat3<real_t> mat_bravais_fcc{ const mat3_t<real_t> mat_bravais_fcc{
0.0, 0.5, 0.0, 0.0, 0.5, 0.0,
0.5, 0.0, 1.0, 0.5, 0.0, 1.0,
0.5, 0.5, 0.0, 0.5, 0.5, 0.0,
}; };
const mat3<real_t> mat_reciprocal_fcc{ const mat3_t<real_t> mat_reciprocal_fcc{
-fourpi, fourpi, twopi, -fourpi, fourpi, twopi,
0.0, 0.0, twopi, 0.0, 0.0, twopi,
fourpi, 0.0, -twopi, fourpi, 0.0, -twopi,
}; };
const mat3<int> mat_invrecip_fcc{ const mat3_t<int> mat_invrecip_fcc{
0, 1, 1, 0, 1, 1,
1, 0, 1, 1, 0, 1,
0, 2, 0, 0, 2, 0,
}; };
const std::vector<vec3<real_t>> normals_fcc{ const std::vector<vec3_t<real_t>> normals_fcc{
{twopi,0.,0.},{-twopi,0.,0.}, {twopi,0.,0.},{-twopi,0.,0.},
{0.,twopi,0.},{0.,-twopi,0.}, {0.,twopi,0.},{0.,-twopi,0.},
{0.,0.,twopi},{0.,0.,-twopi}, {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; }; auto kronecker = []( int i, int j ) -> real_t { return (i==j)? 1.0 : 0.0; };
//! Ewald summation: short-range Green's function //! Ewald summation: short-range Green's function
auto add_greensftide_sr = [&]( mat3<real_t>& D, const vec3<real_t>& d ) -> void { auto add_greensftide_sr = [&]( mat3_t<real_t>& D, const vec3_t<real_t>& d ) -> void {
auto r = d.norm(); auto r = d.norm();
if( r< 1e-14 ) return; // return zero for r=0 if( r< 1e-14 ) return; // return zero for r=0
@ -170,7 +170,7 @@ private:
}; };
//! Ewald summation: long-range Green's function //! Ewald summation: long-range Green's function
auto add_greensftide_lr = [&]( mat3<real_t>& D, const vec3<real_t>& k, const vec3<real_t>& r ) -> void { auto add_greensftide_lr = [&]( mat3_t<real_t>& D, const vec3_t<real_t>& k, const vec3_t<real_t>& r ) -> void {
real_t kmod2 = k.norm_squared(); real_t kmod2 = k.norm_squared();
real_t term = std::exp(-kmod2/(4*alpha2))*std::cos(k.dot(r)) / kmod2 * fft_norm; real_t term = std::exp(-kmod2/(4*alpha2))*std::cos(k.dot(r)) / kmod2 * fft_norm;
for( int mu=0; mu<3; ++mu ){ for( int mu=0; mu<3; ++mu ){
@ -195,22 +195,22 @@ private:
constexpr ptrdiff_t lnumber = 3, knumber = 3; constexpr ptrdiff_t lnumber = 3, knumber = 3;
const int numb = 1; //!< search radius when shifting vectors into FBZ const int numb = 1; //!< search radius when shifting vectors into FBZ
vectk_.assign(D_xx_.memsize(),vec3<real_t>()); vectk_.assign(D_xx_.memsize(),vec3_t<real_t>());
ico_.assign(D_xx_.memsize(),vec3<int>()); ico_.assign(D_xx_.memsize(),vec3_t<int>());
vecitk_.assign(D_xx_.memsize(),vec3<int>()); vecitk_.assign(D_xx_.memsize(),vec3_t<int>());
#pragma omp parallel #pragma omp parallel
{ {
//... temporary to hold values of the dynamical matrix //... temporary to hold values of the dynamical matrix
mat3<real_t> matD(0.0); mat3_t<real_t> matD(0.0);
#pragma omp for #pragma omp for
for( ptrdiff_t i=0; i<nlattice; ++i ){ for( ptrdiff_t i=0; i<nlattice; ++i ){
for( ptrdiff_t j=0; j<nlattice; ++j ){ for( ptrdiff_t j=0; j<nlattice; ++j ){
for( ptrdiff_t k=0; k<nlattice; ++k ){ for( ptrdiff_t k=0; k<nlattice; ++k ){
// compute lattice site vector from (i,j,k) multiplying Bravais base matrix, and wrap back to box // compute lattice site vector from (i,j,k) multiplying Bravais base matrix, and wrap back to box
const vec3<real_t> x_ijk({dx*real_t(i),dx*real_t(j),dx*real_t(k)}); const vec3_t<real_t> x_ijk({dx*real_t(i),dx*real_t(j),dx*real_t(k)});
const vec3<real_t> ar = (mat_bravais * x_ijk).wrap_abs(); const vec3_t<real_t> ar = (mat_bravais * x_ijk).wrap_abs();
//... zero temporary matrix //... zero temporary matrix
matD.zero(); matD.zero();
@ -219,8 +219,8 @@ private:
for( ptrdiff_t ix=-lnumber; ix<=lnumber; ix++ ){ for( ptrdiff_t ix=-lnumber; ix<=lnumber; ix++ ){
for( ptrdiff_t iy=-lnumber; iy<=lnumber; iy++ ){ for( ptrdiff_t iy=-lnumber; iy<=lnumber; iy++ ){
for( ptrdiff_t iz=-lnumber; iz<=lnumber; iz++ ){ for( ptrdiff_t iz=-lnumber; iz<=lnumber; iz++ ){
const vec3<real_t> n_ijk({real_t(ix),real_t(iy),real_t(iz)}); const vec3_t<real_t> n_ijk({real_t(ix),real_t(iy),real_t(iz)});
const vec3<real_t> dr(ar - mat_bravais * n_ijk); const vec3_t<real_t> dr(ar - mat_bravais * n_ijk);
add_greensftide_sr(matD, dr); add_greensftide_sr(matD, dr);
} }
} }
@ -231,8 +231,8 @@ private:
for( ptrdiff_t iy=-knumber; iy<=knumber; iy++ ){ for( ptrdiff_t iy=-knumber; iy<=knumber; iy++ ){
for( ptrdiff_t iz=-knumber; iz<=knumber; iz++ ){ for( ptrdiff_t iz=-knumber; iz<=knumber; iz++ ){
if(std::abs(ix)+std::abs(iy)+std::abs(iz) != 0){ if(std::abs(ix)+std::abs(iy)+std::abs(iz) != 0){
const vec3<real_t> k_ijk({real_t(ix)/nlattice,real_t(iy)/nlattice,real_t(iz)/nlattice}); const vec3_t<real_t> k_ijk({real_t(ix)/nlattice,real_t(iy)/nlattice,real_t(iz)/nlattice});
const vec3<real_t> ak( mat_reciprocal * k_ijk); const vec3_t<real_t> ak( mat_reciprocal * k_ijk);
add_greensftide_lr(matD, ak, ar ); add_greensftide_lr(matD, ak, ar );
} }
@ -278,7 +278,7 @@ private:
std::ofstream ofs2("test_brillouin.txt"); std::ofstream ofs2("test_brillouin.txt");
#endif #endif
using map_t = std::map<vec3<int>,size_t>; using map_t = std::map<vec3_t<int>,size_t>;
map_t iimap; map_t iimap;
//!=== Make temporary copies before resorting to std. Fourier grid ========!// //!=== Make temporary copies before resorting to std. Fourier grid ========!//
@ -312,8 +312,8 @@ private:
#pragma omp parallel #pragma omp parallel
{ {
// thread private matrix representation // thread private matrix representation
mat3<real_t> D; mat3_t<real_t> D;
vec3<real_t> eval, evec1, evec2, evec3; vec3_t<real_t> eval, evec1, evec2, evec3_t;
#pragma omp for #pragma omp for
for( size_t i=0; i<D_xx_.size(0); i++ ) for( size_t i=0; i<D_xx_.size(0); i++ )
@ -322,7 +322,7 @@ private:
{ {
for( size_t k=0; k<D_xx_.size(2); k++ ) for( size_t k=0; k<D_xx_.size(2); k++ )
{ {
vec3<real_t> kv = D_xx_.get_k<real_t>(i,j,k); vec3_t<real_t> kv = D_xx_.get_k<real_t>(i,j,k);
// put matrix elements into actual matrix // put matrix elements into actual matrix
D(0,0) = std::real(temp1.kelem(i,j,k)) / fft_norm12; 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; D(2,2) = std::imag(temp3.kelem(i,j,k)) / fft_norm12;
// compute eigenstructure of matrix // compute eigenstructure of matrix
D.eigen(eval, evec1, evec2, evec3); D.eigen(eval, evec1, evec2, evec3_t);
evec3 /= (twopi*ngrid_); evec3_t /= (twopi*ngrid_);
// now determine to which modes on the regular lattice this contributes // now determine to which modes on the regular lattice this contributes
vec3<real_t> ar = kv / (twopi*ngrid_); vec3_t<real_t> ar = kv / (twopi*ngrid_);
vec3<real_t> a(mat_reciprocal * ar); vec3_t<real_t> a(mat_reciprocal * ar);
// translate the k-vectors into the "candidate" FBZ // translate the k-vectors into the "candidate" FBZ
for( int l1=-numb; l1<=numb; ++l1 ){ for( int l1=-numb; l1<=numb; ++l1 ){
@ -347,9 +347,9 @@ private:
// need both halfs of Fourier space since we use real transforms // need both halfs of Fourier space since we use real transforms
for( int isign=0; isign<=1; ++isign ){ for( int isign=0; isign<=1; ++isign ){
const real_t sign = 2.0*real_t(isign)-1.0; const real_t sign = 2.0*real_t(isign)-1.0;
const vec3<real_t> vshift({real_t(l1),real_t(l2),real_t(l3)}); const vec3_t<real_t> vshift({real_t(l1),real_t(l2),real_t(l3)});
vec3<real_t> vectk = sign * a + mat_reciprocal * vshift; vec3_t<real_t> vectk = sign * a + mat_reciprocal * vshift;
if( check_FBZ( normals, vectk ) ) if( check_FBZ( normals, vectk ) )
{ {
@ -358,11 +358,11 @@ private:
int iz = std::round(vectk.z*(ngrid_)/twopi); int iz = std::round(vectk.z*(ngrid_)/twopi);
#pragma omp critical #pragma omp critical
{iimap.insert( std::pair<vec3<int>,size_t>({ix,iy,iz}, D_xx_.get_idx(i,j,k)) );} {iimap.insert( std::pair<vec3_t<int>,size_t>({ix,iy,iz}, D_xx_.get_idx(i,j,k)) );}
temp1.kelem(i,j,k) = ccomplex_t(eval[2],eval[1]); temp1.kelem(i,j,k) = ccomplex_t(eval[2],eval[1]);
temp2.kelem(i,j,k) = ccomplex_t(eval[0],evec3.x); temp2.kelem(i,j,k) = ccomplex_t(eval[0],evec3_t.x);
temp3.kelem(i,j,k) = ccomplex_t(evec3.y,evec3.z); temp3.kelem(i,j,k) = ccomplex_t(evec3_t.y,evec3_t.z);
} }
}//sign }//sign
} //l3 } //l3
@ -389,24 +389,24 @@ private:
int ii = (int(i)>nlattice/2)? int(i)-nlattice : int(i); int ii = (int(i)>nlattice/2)? int(i)-nlattice : int(i);
int jj = (int(j)>nlattice/2)? int(j)-nlattice : int(j); int jj = (int(j)>nlattice/2)? int(j)-nlattice : int(j);
int kk = (int(k)>nlattice/2)? int(k)-nlattice : int(k); int kk = (int(k)>nlattice/2)? int(k)-nlattice : int(k);
vec3<real_t> kv({real_t(ii),real_t(jj),real_t(kk)}); vec3_t<real_t> kv({real_t(ii),real_t(jj),real_t(kk)});
auto align_with_k = [&]( const vec3<real_t>& v ) -> vec3<real_t>{ auto align_with_k = [&]( const vec3_t<real_t>& v ) -> vec3_t<real_t>{
return v*((v.dot(kv)<0.0)?-1.0:1.0); return v*((v.dot(kv)<0.0)?-1.0:1.0);
}; };
vec3<real_t> v, l; vec3_t<real_t> v, l;
map_t::iterator it; map_t::iterator it;
if( !is_in(i,j,k,mat_invrecip) ){ if( !is_in(i,j,k,mat_invrecip) ){
auto average_lv = [&]( const auto& t1, const auto& t2, const auto& t3, vec3<real_t>& v, vec3<real_t>& l ) { auto average_lv = [&]( const auto& t1, const auto& t2, const auto& t3, vec3_t<real_t>& v, vec3_t<real_t>& l ) {
v = 0.0; l = 0.0; v = 0.0; l = 0.0;
int count(0); int count(0);
auto add_lv = [&]( auto it ) -> void { auto add_lv = [&]( auto it ) -> void {
auto q = it->second;++count; auto q = it->second;++count;
l += vec3<real_t>({std::real(t1.kelem(q)),std::imag(t1.kelem(q)),std::real(t2.kelem(q))}); l += vec3_t<real_t>({std::real(t1.kelem(q)),std::imag(t1.kelem(q)),std::real(t2.kelem(q))});
v += align_with_k(vec3<real_t>({std::imag(t2.kelem(q)),std::real(t3.kelem(q)),std::imag(t3.kelem(q))})); v += align_with_k(vec3_t<real_t>({std::imag(t2.kelem(q)),std::real(t3.kelem(q)),std::imag(t3.kelem(q))}));
}; };
map_t::iterator it; map_t::iterator it;
if( (it = iimap.find({ii-1,jj,kk}))!=iimap.end() ){ add_lv(it); } if( (it = iimap.find({ii-1,jj,kk}))!=iimap.end() ){ add_lv(it); }
@ -423,8 +423,8 @@ private:
}else{ }else{
if( (it = iimap.find({ii,jj,kk}))!=iimap.end() ){ if( (it = iimap.find({ii,jj,kk}))!=iimap.end() ){
auto q = it->second; auto q = it->second;
l = vec3<real_t>({std::real(temp1.kelem(q)),std::imag(temp1.kelem(q)),std::real(temp2.kelem(q))}); l = vec3_t<real_t>({std::real(temp1.kelem(q)),std::imag(temp1.kelem(q)),std::real(temp2.kelem(q))});
v = align_with_k(vec3<real_t>({std::imag(temp2.kelem(q)),std::real(temp3.kelem(q)),std::imag(temp3.kelem(q))})); v = align_with_k(vec3_t<real_t>({std::imag(temp2.kelem(q)),std::real(temp3.kelem(q)),std::imag(temp3.kelem(q))}));
} }
} }
D_xx_.kelem(i,j,k) = l[0]; D_xx_.kelem(i,j,k) = l[0];
@ -443,13 +443,13 @@ private:
for( size_t j=0; j<D_xx_.size(1); j++ ){ for( size_t j=0; j<D_xx_.size(1); j++ ){
for( size_t k=0; k<D_xx_.size(2); k++ ) for( size_t k=0; k<D_xx_.size(2); k++ )
{ {
vec3<real_t> kv = D_xx_.get_k<real_t>(i,j,k); vec3_t<real_t> kv = D_xx_.get_k<real_t>(i,j,k);
double mu1 = std::real(D_xx_.kelem(i,j,k)); double mu1 = std::real(D_xx_.kelem(i,j,k));
// double mu2 = std::real(D_xy_.kelem(i,j,k)); // double mu2 = std::real(D_xy_.kelem(i,j,k));
// double mu3 = std::real(D_xz_.kelem(i,j,k)); // double mu3 = std::real(D_xz_.kelem(i,j,k));
vec3<real_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))}); vec3_t<real_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(); 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 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); real_t st = std::sin(ktheta), ct = std::cos(ktheta), sp = std::sin(kphi), cp = std::cos(kphi);
vec3<real_t> e_r( st*cp, st*sp, ct), e_theta( ct*cp, ct*sp, -st), e_phi( -sp, cp, 0.0 ); vec3_t<real_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 // re-normalise to that longitudinal amplitude is exact
double renorm = evec1.dot( e_r ); if( renorm < 0.01 ) renorm = 1.0; double renorm = evec1.dot( e_r ); if( renorm < 0.01 ) renorm = 1.0;

56
include/vec3.hh
View file

@ -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 a code to generate initial conditions for cosmological simulations
CHANGELOG (only majors, for details see repo): CHANGELOG (only majors, for details see repo):
@ -9,7 +9,7 @@
//! implements a simple class of 3-vectors of arbitrary scalar type //! implements a simple class of 3-vectors of arbitrary scalar type
template< typename T > template< typename T >
class vec3{ class vec3_t{
private: private:
//! holds the data //! holds the data
std::array<T,3> data_; std::array<T,3> data_;
@ -19,27 +19,27 @@ public:
T &x,&y,&z; T &x,&y,&z;
//! empty constructor //! empty constructor
vec3() vec3_t()
: x(data_[0]),y(data_[1]),z(data_[2]){} : x(data_[0]),y(data_[1]),z(data_[2]){}
//! copy constructor //! copy constructor
vec3( const vec3<T> &v) vec3_t( const vec3_t<T> &v)
: data_(v.data_), x(data_[0]),y(data_[1]),z(data_[2]){} : 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 //! copy constructor for non-const reference, needed to avoid variadic template being called for non-const reference
vec3( vec3<T>& v) vec3_t( vec3_t<T>& v)
: data_(v.data_), x(data_[0]),y(data_[1]),z(data_[2]){} : data_(v.data_), x(data_[0]),y(data_[1]),z(data_[2]){}
//! move constructor //! move constructor
vec3( vec3<T> &&v) vec3_t( vec3_t<T> &&v)
: data_(std::move(v.data_)), x(data_[0]), y(data_[1]), z(data_[2]){} : 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<typename ...E> template<typename ...E>
vec3(E&&...e) vec3_t(E&&...e)
: data_{{std::forward<E>(e)...}}, x{data_[0]}, y{data_[1]}, z{data_[2]} : data_{{std::forward<E>(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]){} // : data_{{a,b,c}}, x(data_[0]), y(data_[1]), z(data_[2]){}
//! bracket index access to vector components //! bracket index access to vector components
@ -49,37 +49,37 @@ public:
const T &operator[](size_t i) const noexcept { return data_[i]; } const T &operator[](size_t i) const noexcept { return data_[i]; }
// assignment operator // assignment operator
vec3<T>& operator=( const vec3<T>& v ) noexcept { data_=v.data_; return *this; } vec3_t<T>& operator=( const vec3_t<T>& v ) noexcept { data_=v.data_; return *this; }
//! implementation of summation of vec3 //! 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}); } vec3_t<T> operator+( const vec3_t<T>& v ) const noexcept{ return vec3_t<T>({x+v.x,y+v.y,z+v.z}); }
//! implementation of difference of vec3 //! 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}); } vec3_t<T> operator-( const vec3_t<T>& v ) const noexcept{ return vec3_t<T>({x-v.x,y-v.y,z-v.z}); }
//! implementation of unary negative //! implementation of unary negative
vec3<T> operator-() const noexcept{ return vec3<T>({-x,-y,-z}); } vec3_t<T> operator-() const noexcept{ return vec3_t<T>({-x,-y,-z}); }
//! implementation of scalar multiplication //! implementation of scalar multiplication
vec3<T> operator*( T s ) const noexcept{ return vec3<T>({x*s,y*s,z*s}); } vec3_t<T> operator*( T s ) const noexcept{ return vec3_t<T>({x*s,y*s,z*s}); }
//! implementation of scalar division //! implementation of scalar division
vec3<T> operator/( T s ) const noexcept{ return vec3<T>({x/s,y/s,z/s}); } vec3_t<T> operator/( T s ) const noexcept{ return vec3_t<T>({x/s,y/s,z/s}); }
//! implementation of += operator //! implementation of += operator
vec3<T>& operator+=( const vec3<T>& v ) noexcept{ x+=v.x; y+=v.y; z+=v.z; return *this; } vec3_t<T>& operator+=( const vec3_t<T>& v ) noexcept{ x+=v.x; y+=v.y; z+=v.z; return *this; }
//! implementation of -= operator //! implementation of -= operator
vec3<T>& operator-=( const vec3<T>& v ) noexcept{ x-=v.x; y-=v.y; z-=v.z; return *this; } vec3_t<T>& operator-=( const vec3_t<T>& v ) noexcept{ x-=v.x; y-=v.y; z-=v.z; return *this; }
//! multiply with scalar //! multiply with scalar
vec3<T>& operator*=( T s ) noexcept{ x*=s; y*=s; z*=s; return *this; } vec3_t<T>& operator*=( T s ) noexcept{ x*=s; y*=s; z*=s; return *this; }
//! divide by scalar //! divide by scalar
vec3<T>& operator/=( T s ) noexcept{ x/=s; y/=s; z/=s; return *this; } vec3_t<T>& operator/=( T s ) noexcept{ x/=s; y/=s; z/=s; return *this; }
//! compute dot product with another vector //! compute dot product with another vector
T dot(const vec3<T> &a) const noexcept T dot(const vec3_t<T> &a) const noexcept
{ {
return data_[0] * a.data_[0] + data_[1] * a.data_[1] + data_[2] * a.data_[2]; 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() ); } T norm(void) const noexcept { return std::sqrt( this->norm_squared() ); }
//! wrap absolute vector to box of size p //! wrap absolute vector to box of size p
vec3<T>& wrap_abs( T p = 1.0 ) noexcept{ vec3_t<T>& wrap_abs( T p = 1.0 ) noexcept{
for( auto& x : data_ ) x = std::fmod( 2*p + x, p ); for( auto& x : data_ ) x = std::fmod( 2*p + x, p );
return *this; return *this;
} }
//! wrap relative vector to box of size p //! wrap relative vector to box of size p
vec3<T>& wrap_rel( T p = 1.0 ) noexcept{ vec3_t<T>& wrap_rel( T p = 1.0 ) noexcept{
for( auto& x : data_ ) x = (x<-p/2)? x+p : (x>=p/2)? x-p : x; for( auto& x : data_ ) x = (x<-p/2)? x+p : (x>=p/2)? x-p : x;
return *this; return *this;
} }
//! ordering, allows 3d sorting of vec3s //! ordering, allows 3d sorting of vec3_ts
bool operator<( const vec3<T>& o ) const noexcept{ bool operator<( const vec3_t<T>& o ) const noexcept{
if( x!=o.x ) return x<o.x?true:false; if( x!=o.x ) return x<o.x?true:false;
if( y!=o.y ) return y<o.y?true:false; if( y!=o.y ) return y<o.y?true:false;
if( z!=o.z ) return z<o.z?true:false; if( z!=o.z ) return z<o.z?true:false;
@ -113,6 +113,6 @@ public:
//! multiplication with scalar //! multiplication with scalar
template<typename T> template<typename T>
vec3<T> operator*( T s, const vec3<T>& v ){ vec3_t<T> operator*( T s, const vec3_t<T>& v ){
return vec3<T>({v.x*s,v.y*s,v.z*s}); return vec3_t<T>({v.x*s,v.y*s,v.z*s});
} }

2
src/grid_fft.cc
View file

@ -860,7 +860,7 @@ void Grid_FFT<data_t,bdistributed>::Compute_PowerSpectrum(std::vector<double> &b
for (size_t iy = 0; iy < size(1); iy++) for (size_t iy = 0; iy < size(1); iy++)
for (size_t iz = 0; iz < size(2); iz++) for (size_t iz = 0; iz < size(2); iz++)
{ {
vec3<double> k3 = get_k<double>(ix, iy, iz); vec3_t<double> k3 = get_k<double>(ix, iy, iz);
double k = k3.norm(); double k = k3.norm();
int idx2 = k / dk; //int((1.0f / dklog * std::log10(k / kmin))); int idx2 = k / dk; //int((1.0f / dklog * std::log10(k / kmin)));
auto z = this->kelem(ix, iy, iz); auto z = this->kelem(ix, iy, iz);

6
src/ic_generator.cc
View file

@ -320,7 +320,7 @@ int Run( ConfigFile& the_config )
if (bAddExternalTides) if (bAddExternalTides)
{ {
phi2.assign_function_of_grids_kdep([&](vec3<real_t> kvec, ccomplex_t pphi, ccomplex_t pphi2) { phi2.assign_function_of_grids_kdep([&](vec3_t<real_t> kvec, ccomplex_t pphi, ccomplex_t pphi2) {
// sign in front of f_aniso is reversed since phi1 = -phi // 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; 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) ); + 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 ){ if( bDoBaryons ){
vec3<real_t> kvec = phi.get_k<real_t>(i,j,k); vec3_t<real_t> kvec = phi.get_k<real_t>(i,j,k);
real_t k2 = kvec.norm_squared(), kmod = std::sqrt(k2); real_t k2 = kvec.norm_squared(), kmod = std::sqrt(k2);
double ampldiff = ((this_species == cosmo_species::dm)? the_cosmo_calc->GetAmplitude(kmod, cdm0) : double ampldiff = ((this_species == cosmo_species::dm)? the_cosmo_calc->GetAmplitude(kmod, cdm0) :
(this_species == cosmo_species::baryon)? the_cosmo_calc->GetAmplitude(kmod, baryon0) : (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)) ); + 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 ){ if( bDoBaryons ){
vec3<real_t> kvec = phi.get_k<real_t>(i,j,k); vec3_t<real_t> kvec = phi.get_k<real_t>(i,j,k);
real_t k2 = kvec.norm_squared(), kmod = std::sqrt(k2); real_t k2 = kvec.norm_squared(), kmod = std::sqrt(k2);
double ampldiff = ((this_species == cosmo_species::dm)? the_cosmo_calc->GetAmplitude(kmod, vcdm0) : double ampldiff = ((this_species == cosmo_species::dm)? the_cosmo_calc->GetAmplitude(kmod, vcdm0) :
(this_species == cosmo_species::baryon)? the_cosmo_calc->GetAmplitude(kmod, vbaryon0) : (this_species == cosmo_species::baryon)? the_cosmo_calc->GetAmplitude(kmod, vbaryon0) :

3
src/plugins/transfer_CLASS.cc
View file

@ -223,13 +223,14 @@ public:
gsl_spline *splineT = nullptr; gsl_spline *splineT = nullptr;
gsl_interp_accel *accT = nullptr; gsl_interp_accel *accT = nullptr;
switch(type){ switch(type){
// values at ztarget:
case total: splineT = gsl_sp_dtot_; accT = gsl_ia_dtot_; break; case total: splineT = gsl_sp_dtot_; accT = gsl_ia_dtot_; break;
case cdm: splineT = gsl_sp_dc_; accT = gsl_ia_dc_; break; case cdm: splineT = gsl_sp_dc_; accT = gsl_ia_dc_; break;
case baryon: splineT = gsl_sp_db_; accT = gsl_ia_db_; break; case baryon: splineT = gsl_sp_db_; accT = gsl_ia_db_; break;
case vtotal: splineT = gsl_sp_ttot_; accT = gsl_ia_ttot_; break; case vtotal: splineT = gsl_sp_ttot_; accT = gsl_ia_ttot_; break;
case vcdm: splineT = gsl_sp_tc_; accT = gsl_ia_tc_; break; case vcdm: splineT = gsl_sp_tc_; accT = gsl_ia_tc_; break;
case vbaryon: splineT = gsl_sp_tb_; accT = gsl_ia_tb_; 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 total0: splineT = gsl_sp_dtot0_;accT = gsl_ia_dtot0_;break;
case cdm0: splineT = gsl_sp_dc0_; accT = gsl_ia_dc0_; break; case cdm0: splineT = gsl_sp_dc0_; accT = gsl_ia_dc0_; break;
case baryon0: splineT = gsl_sp_db0_; accT = gsl_ia_db0_; break; case baryon0: splineT = gsl_sp_db0_; accT = gsl_ia_db0_; break;