mirror of
https://github.com/cosmo-sims/monofonIC.git
synced 2024-09-19 17:03:45 +02:00
fixed compilation problems when changing code precision
This commit is contained in:
Oliver Hahn
parent
3bead20876
commit
326bfb7b9e
10 changed files with 106 additions and 110 deletions
|
|
@ -84,9 +84,9 @@ public:
|
|||
|
||||
void reset()
|
||||
{
|
||||
if (data_ != nullptr) { fftw_free(data_); data_ = nullptr; }
|
||||
if (plan_ != nullptr) { fftw_destroy_plan(plan_); plan_ = nullptr; }
|
||||
if (iplan_ != nullptr) { fftw_destroy_plan(iplan_); iplan_ = nullptr; }
|
||||
if (data_ != nullptr) { FFTW_API(free)(data_); data_ = nullptr; }
|
||||
if (plan_ != nullptr) { FFTW_API(destroy_plan)(plan_); plan_ = nullptr; }
|
||||
if (iplan_ != nullptr) { FFTW_API(destroy_plan)(iplan_); iplan_ = nullptr; }
|
||||
ballocated_ = false;
|
||||
}
|
||||
|
||||
|
|
@ -305,9 +305,9 @@ public:
|
|||
data_t get_cic( const vec3_t<real_t>& v ) const noexcept
|
||||
{
|
||||
// warning! this doesn't work with MPI
|
||||
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.z/length_[2]+1.0,1.0)*n_[2] });
|
||||
vec3_t<real_t> x({real_t(std::fmod(v.x/length_[0]+1.0,1.0)*n_[0]),
|
||||
real_t(std::fmod(v.y/length_[1]+1.0,1.0)*n_[1]),
|
||||
real_t(std::fmod(v.z/length_[2]+1.0,1.0)*n_[2]) });
|
||||
size_t ix = static_cast<size_t>(x.x);
|
||||
size_t iy = static_cast<size_t>(x.y);
|
||||
size_t iz = static_cast<size_t>(x.z);
|
||||
|
|
|
|||
|
|
@ -27,6 +27,7 @@ template<typename T>
|
|||
class mat3_t{
|
||||
protected:
|
||||
std::array<T,9> data_;
|
||||
std::array<double,9> data_double_;
|
||||
gsl_matrix_view m_;
|
||||
gsl_vector *eval_;
|
||||
gsl_matrix *evec_;
|
||||
|
|
@ -37,12 +38,20 @@ protected:
|
|||
// allocate memory for GSL operations if we haven't done so yet
|
||||
if( !bdid_alloc_gsl_ )
|
||||
{
|
||||
m_ = gsl_matrix_view_array (&data_[0], 3, 3);
|
||||
if( typeid(T)!=typeid(double) ){
|
||||
m_ = gsl_matrix_view_array ( &data_double_[0], 3, 3);
|
||||
}else{
|
||||
m_ = gsl_matrix_view_array ( (double*)(&data_[0]), 3, 3); // this should only ever be called for T==double so cast is to avoid constexpr ifs from C++17
|
||||
}
|
||||
eval_ = gsl_vector_alloc (3);
|
||||
evec_ = gsl_matrix_alloc (3, 3);
|
||||
wsp_ = gsl_eigen_symmv_alloc (3);
|
||||
bdid_alloc_gsl_ = true;
|
||||
}
|
||||
|
||||
if( typeid(T)!=typeid(double) ){
|
||||
for( int i=0; i<9; ++i ) data_double_[i] = double(data_[i]);
|
||||
}
|
||||
}
|
||||
|
||||
void free_gsl(){
|
||||
|
|
|
|||
|
|
@ -39,19 +39,19 @@ namespace particle
|
|||
const std::vector<std::vector<vec3_t<real_t>>> lattice_shifts =
|
||||
{
|
||||
// first shift must always be zero! (otherwise set_positions and set_velocities break)
|
||||
/* SC : */ {{0.0, 0.0, 0.0}},
|
||||
/* BCC: */ {{0.0, 0.0, 0.0}, {0.5, 0.5, 0.5}},
|
||||
/* FCC: */ {{0.0, 0.0, 0.0}, {0.0, 0.5, 0.5}, {0.5, 0.0, 0.5}, {0.5, 0.5, 0.0}},
|
||||
/* 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}},
|
||||
/* SC : */ {{real_t{0.0}, real_t{0.0}, real_t{0.0}}},
|
||||
/* BCC: */ {{real_t{0.0}, real_t{0.0}, real_t{0.0}}, {real_t{0.5}, real_t{0.5}, real_t{0.5}}},
|
||||
/* FCC: */ {{real_t{0.0}, real_t{0.0}, real_t{0.0}}, {real_t{0.0}, real_t{0.5}, real_t{0.5}}, {real_t{0.5}, real_t{0.0}, real_t{0.5}}, {real_t{0.5}, real_t{0.5}, real_t{0.0}}},
|
||||
/* RSC: */ {{real_t{0.0}, real_t{0.0}, real_t{0.0}}, {real_t{0.0}, real_t{0.0}, real_t{0.5}}, {real_t{0.0}, real_t{0.5}, real_t{0.0}}, {real_t{0.0}, real_t{0.5}, real_t{0.5}}, {real_t{0.5}, real_t{0.0}, real_t{0.0}}, {real_t{0.5}, real_t{0.0}, real_t{0.5}}, {real_t{0.5}, real_t{0.5}, real_t{0.0}}, {real_t{0.5}, real_t{0.5}, real_t{0.5}}},
|
||||
};
|
||||
|
||||
const std::vector<vec3_t<real_t>> 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?!?
|
||||
/* FCC: */ {0.5, 0.5, 0.5}, // this corresponds to NaCl lattice
|
||||
// /* FCC: */ {0.25, 0.25, 0.25}, // this corresponds to Zincblende/GaAs lattice
|
||||
/* RSC: */ {0.25, 0.25, 0.25},
|
||||
/* SC : */ {real_t{0.5}, real_t{0.5}, real_t{0.5}}, // this corresponds to CsCl lattice
|
||||
/* BCC: */ {real_t{0.5}, real_t{0.5}, real_t{0.0}}, // is there a diatomic lattice with BCC base?!?
|
||||
/* FCC: */ {real_t{0.5}, real_t{0.5}, real_t{0.5}}, // this corresponds to NaCl lattice
|
||||
// /* FCC: */ {real_t{0.25}, real_t{0.25}, real_t{0.25}}, // this corresponds to Zincblende/GaAs lattice
|
||||
/* RSC: */ {real_t{0.25}, real_t{0.25}, real_t{0.25}},
|
||||
};
|
||||
|
||||
template <typename field_t>
|
||||
|
|
@ -283,7 +283,7 @@ namespace particle
|
|||
{
|
||||
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] + (is_second_lattice ? second_lattice_shift[lattice_type] : vec3_t<real_t>{0., 0., 0.}));
|
||||
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_t<real_t>{real_t(0.), real_t(0.), real_t(0.)}));
|
||||
if (b64reals)
|
||||
{
|
||||
particles_.set_pos64(ipcount++, idim, pos[idim] * lunit + field.relem(i, j, k));
|
||||
|
|
|
|||
|
|
@ -76,14 +76,14 @@ private:
|
|||
//! === vectors, reciprocals and normals for the SC lattice ===
|
||||
const int charge_fac_sc = 1;
|
||||
const mat3_t<real_t> mat_bravais_sc{
|
||||
1.0, 0.0, 0.0,
|
||||
0.0, 1.0, 0.0,
|
||||
0.0, 0.0, 1.0,
|
||||
real_t{1.0}, real_t{0.0}, real_t{0.0},
|
||||
real_t{0.0}, real_t{1.0}, real_t{0.0},
|
||||
real_t{0.0}, real_t{0.0}, real_t{1.0},
|
||||
};
|
||||
const mat3_t<real_t> mat_reciprocal_sc{
|
||||
twopi, 0.0, 0.0,
|
||||
0.0, twopi, 0.0,
|
||||
0.0, 0.0, twopi,
|
||||
twopi, real_t{0.0}, real_t{0.0},
|
||||
real_t{0.0}, twopi, real_t{0.0},
|
||||
real_t{0.0}, real_t{0.0}, twopi,
|
||||
};
|
||||
const mat3_t<int> mat_invrecip_sc{
|
||||
2, 0, 0,
|
||||
|
|
@ -91,22 +91,22 @@ private:
|
|||
0, 0, 2,
|
||||
};
|
||||
const std::vector<vec3_t<real_t>> normals_sc{
|
||||
{pi,0.,0.},{-pi,0.,0.},
|
||||
{0.,pi,0.},{0.,-pi,0.},
|
||||
{0.,0.,pi},{0.,0.,-pi},
|
||||
{pi,real_t{0.},real_t{0.}},{-pi,real_t{0.},real_t{0.}},
|
||||
{real_t{0.},pi,real_t{0.}},{real_t{0.},-pi,real_t{0.}},
|
||||
{real_t{0.},real_t{0.},pi},{real_t{0.},real_t{0.},-pi},
|
||||
};
|
||||
|
||||
|
||||
//! === vectors, reciprocals and normals for the BCC lattice ===
|
||||
const int charge_fac_bcc = 2;
|
||||
const mat3_t<real_t> mat_bravais_bcc{
|
||||
1.0, 0.0, 0.5,
|
||||
0.0, 1.0, 0.5,
|
||||
0.0, 0.0, 0.5,
|
||||
real_t{1.0}, real_t{0.0}, real_t{0.5},
|
||||
real_t{0.0}, real_t{1.0}, real_t{0.5},
|
||||
real_t{0.0}, real_t{0.0}, real_t{0.5},
|
||||
};
|
||||
const mat3_t<real_t> mat_reciprocal_bcc{
|
||||
twopi, 0.0, 0.0,
|
||||
0.0, twopi, 0.0,
|
||||
twopi, real_t{0.0}, real_t{0.0},
|
||||
real_t{0.0}, twopi, real_t{0.0},
|
||||
-twopi, -twopi, fourpi,
|
||||
};
|
||||
const mat3_t<int> mat_invrecip_bcc{
|
||||
|
|
@ -115,23 +115,23 @@ private:
|
|||
1, 1, 1,
|
||||
};
|
||||
const std::vector<vec3_t<real_t>> 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.}
|
||||
{real_t{0.0},pi,pi},{real_t{0.0},-pi,pi},{real_t{0.0},pi,-pi},{real_t{0.0},-pi,-pi},
|
||||
{pi,real_t{0.0},pi},{-pi,real_t{0.0},pi},{pi,real_t{0.0},-pi},{-pi,real_t{0.0},-pi},
|
||||
{pi,pi,real_t{0.0}},{-pi,pi,real_t{0.0}},{pi,-pi,real_t{0.0}},{-pi,-pi,real_t{0.0}}
|
||||
};
|
||||
|
||||
|
||||
//! === vectors, reciprocals and normals for the FCC lattice ===
|
||||
const int charge_fac_fcc = 4;
|
||||
const mat3_t<real_t> mat_bravais_fcc{
|
||||
0.0, 0.5, 0.0,
|
||||
0.5, 0.0, 1.0,
|
||||
0.5, 0.5, 0.0,
|
||||
real_t{0.0}, real_t{0.5}, real_t{0.0},
|
||||
real_t{0.5}, real_t{0.0}, real_t{1.0},
|
||||
real_t{0.5}, real_t{0.5}, real_t{0.0},
|
||||
};
|
||||
const mat3_t<real_t> mat_reciprocal_fcc{
|
||||
-fourpi, fourpi, twopi,
|
||||
0.0, 0.0, twopi,
|
||||
fourpi, 0.0, -twopi,
|
||||
real_t{0.0}, real_t{0.0}, twopi,
|
||||
fourpi, real_t{0.0}, -twopi,
|
||||
};
|
||||
const mat3_t<int> mat_invrecip_fcc{
|
||||
0, 1, 1,
|
||||
|
|
@ -139,9 +139,9 @@ private:
|
|||
0, 2, 0,
|
||||
};
|
||||
const std::vector<vec3_t<real_t>> normals_fcc{
|
||||
{twopi,0.,0.},{-twopi,0.,0.},
|
||||
{0.,twopi,0.},{0.,-twopi,0.},
|
||||
{0.,0.,twopi},{0.,0.,-twopi},
|
||||
{twopi,real_t{0.0},real_t{0.0}},{-twopi,real_t{0.0},real_t{0.0}},
|
||||
{real_t{0.0},twopi,real_t{0.0}},{real_t{0.0},-twopi,real_t{0.0}},
|
||||
{real_t{0.0},real_t{0.0},twopi},{real_t{0.0},real_t{0.0},-twopi},
|
||||
{+pi,+pi,+pi},{+pi,+pi,-pi},
|
||||
{+pi,-pi,+pi},{+pi,-pi,-pi},
|
||||
{-pi,+pi,+pi},{-pi,+pi,-pi},
|
||||
|
|
@ -223,7 +223,7 @@ private:
|
|||
#pragma omp parallel
|
||||
{
|
||||
//... temporary to hold values of the dynamical matrix
|
||||
mat3_t<real_t> matD(0.0);
|
||||
mat3_t<real_t> matD(real_t(0.0));
|
||||
|
||||
#pragma omp for
|
||||
for( ptrdiff_t i=0; i<nlattice; ++i ){
|
||||
|
|
@ -421,7 +421,7 @@ private:
|
|||
|
||||
if( !is_in(i,j,k,mat_invrecip) ){
|
||||
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 = real_t(0.0); l = real_t(0.0);
|
||||
int count(0);
|
||||
|
||||
auto add_lv = [&]( auto it ) -> void {
|
||||
|
|
@ -476,9 +476,9 @@ private:
|
|||
// ///////////////////////////////////
|
||||
// // project onto spherical coordinate vectors
|
||||
|
||||
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) : real_t(0.0), ktheta = kr>0.0? std::acos( kv.z / kr ): real_t(0.0);
|
||||
real_t st = std::sin(ktheta), ct = std::cos(ktheta), sp = std::sin(kphi), cp = std::cos(kphi);
|
||||
vec3_t<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, real_t(0.0) );
|
||||
|
||||
// re-normalise to that longitudinal amplitude is exact
|
||||
double renorm = evec1.dot( e_r ); if( renorm < 0.01 ) renorm = 1.0;
|
||||
|
|
|
|||
|
|
@ -29,7 +29,7 @@ void Grid_FFT<data_t, bdistributed>::allocate(void)
|
|||
music::dlog.Print("[FFT] Setting up a shared memory field %lux%lux%lu\n", n_[0], n_[1], n_[2]);
|
||||
if (typeid(data_t) == typeid(real_t))
|
||||
{
|
||||
data_ = reinterpret_cast<data_t *>(fftw_malloc(ntot_ * sizeof(real_t)));
|
||||
data_ = reinterpret_cast<data_t *>(FFTW_API(malloc)(ntot_ * sizeof(real_t)));
|
||||
cdata_ = reinterpret_cast<ccomplex_t *>(data_);
|
||||
|
||||
plan_ = FFTW_API(plan_dft_r2c_3d)(n_[0], n_[1], n_[2], (real_t *)data_, (complex_t *)data_, FFTW_RUNMODE);
|
||||
|
|
@ -37,7 +37,7 @@ void Grid_FFT<data_t, bdistributed>::allocate(void)
|
|||
}
|
||||
else if (typeid(data_t) == typeid(ccomplex_t))
|
||||
{
|
||||
data_ = reinterpret_cast<data_t *>(fftw_malloc(ntot_ * sizeof(ccomplex_t)));
|
||||
data_ = reinterpret_cast<data_t *>(FFTW_API(malloc)(ntot_ * sizeof(ccomplex_t)));
|
||||
cdata_ = reinterpret_cast<ccomplex_t *>(data_);
|
||||
|
||||
plan_ = FFTW_API(plan_dft_3d)(n_[0], n_[1], n_[2], (complex_t *)data_, (complex_t *)data_, FFTW_FORWARD, FFTW_RUNMODE);
|
||||
|
|
@ -102,7 +102,7 @@ void Grid_FFT<data_t, bdistributed>::allocate(void)
|
|||
cmplxsz = FFTW_API(mpi_local_size_3d_transposed)(n_[0], n_[1], n_[2] / 2 + 1, MPI_COMM_WORLD,
|
||||
&local_0_size_, &local_0_start_, &local_1_size_, &local_1_start_);
|
||||
ntot_ = 2 * cmplxsz;
|
||||
data_ = (data_t *)fftw_malloc(ntot_ * sizeof(real_t));
|
||||
data_ = (data_t *)FFTW_API(malloc)(ntot_ * sizeof(real_t));
|
||||
cdata_ = reinterpret_cast<ccomplex_t *>(data_);
|
||||
plan_ = FFTW_API(mpi_plan_dft_r2c_3d)(n_[0], n_[1], n_[2], (real_t *)data_, (complex_t *)data_,
|
||||
MPI_COMM_WORLD, FFTW_RUNMODE | FFTW_MPI_TRANSPOSED_OUT);
|
||||
|
|
@ -114,7 +114,7 @@ void Grid_FFT<data_t, bdistributed>::allocate(void)
|
|||
cmplxsz = FFTW_API(mpi_local_size_3d_transposed)(n_[0], n_[1], n_[2], MPI_COMM_WORLD,
|
||||
&local_0_size_, &local_0_start_, &local_1_size_, &local_1_start_);
|
||||
ntot_ = cmplxsz;
|
||||
data_ = (data_t *)fftw_malloc(ntot_ * sizeof(ccomplex_t));
|
||||
data_ = (data_t *)FFTW_API(malloc)(ntot_ * sizeof(ccomplex_t));
|
||||
cdata_ = reinterpret_cast<ccomplex_t *>(data_);
|
||||
plan_ = FFTW_API(mpi_plan_dft_3d)(n_[0], n_[1], n_[2], (complex_t *)data_, (complex_t *)data_,
|
||||
MPI_COMM_WORLD, FFTW_FORWARD, FFTW_RUNMODE | FFTW_MPI_TRANSPOSED_OUT);
|
||||
|
|
|
|||
|
|
@ -140,9 +140,9 @@ int run( config_file& the_config )
|
|||
|
||||
// Anisotropy parameters for beyond box tidal field
|
||||
const std::array<real_t,3> lss_aniso_lambda = {
|
||||
the_config.get_value_safe<double>("cosmology", "LSS_aniso_lx", 0.0),
|
||||
the_config.get_value_safe<double>("cosmology", "LSS_aniso_ly", 0.0),
|
||||
the_config.get_value_safe<double>("cosmology", "LSS_aniso_lz", 0.0),
|
||||
real_t(the_config.get_value_safe<double>("cosmology", "LSS_aniso_lx", 0.0)),
|
||||
real_t(the_config.get_value_safe<double>("cosmology", "LSS_aniso_ly", 0.0)),
|
||||
real_t(the_config.get_value_safe<double>("cosmology", "LSS_aniso_lz", 0.0)),
|
||||
};
|
||||
|
||||
const real_t lss_aniso_sum_lambda = lss_aniso_lambda[0]+lss_aniso_lambda[1]+lss_aniso_lambda[2];
|
||||
|
|
@ -164,24 +164,24 @@ int run( config_file& the_config )
|
|||
const real_t Dplus0 = the_cosmo_calc->get_growth_factor(astart);
|
||||
const real_t vfac = the_cosmo_calc->get_vfact(astart);
|
||||
|
||||
const double g1 = -Dplus0;
|
||||
const double g2 = ((LPTorder>1)? -3.0/7.0*Dplus0*Dplus0 : 0.0);
|
||||
const double g3 = ((LPTorder>2)? 1.0/3.0*Dplus0*Dplus0*Dplus0 : 0.0);
|
||||
const double g3c = ((LPTorder>2)? 1.0/7.0*Dplus0*Dplus0*Dplus0 : 0.0);
|
||||
const real_t g1 = -Dplus0;
|
||||
const real_t g2 = ((LPTorder>1)? -3.0/7.0*Dplus0*Dplus0 : 0.0);
|
||||
const real_t g3 = ((LPTorder>2)? 1.0/3.0*Dplus0*Dplus0*Dplus0 : 0.0);
|
||||
const real_t g3c = ((LPTorder>2)? 1.0/7.0*Dplus0*Dplus0*Dplus0 : 0.0);
|
||||
|
||||
// vfac = d log D+ / dt
|
||||
// d(D+^2)/dt = 2*D+ * d D+/dt = 2 * D+^2 * vfac
|
||||
// d(D+^3)/dt = 3*D+^2* d D+/dt = 3 * D+^3 * vfac
|
||||
const double vfac1 = vfac;
|
||||
const double vfac2 = 2*vfac;
|
||||
const double vfac3 = 3*vfac;
|
||||
const real_t vfac1 = vfac;
|
||||
const real_t vfac2 = 2*vfac;
|
||||
const real_t vfac3 = 3*vfac;
|
||||
|
||||
// anisotropic velocity growth factor for external tides
|
||||
// cf. eq. (5) of Stuecker et al. 2020 (https://arxiv.org/abs/2003.06427)
|
||||
const std::array<real_t,3> lss_aniso_alpha = {
|
||||
1.0 - Dplus0 * lss_aniso_lambda[0],
|
||||
1.0 - Dplus0 * lss_aniso_lambda[1],
|
||||
1.0 - Dplus0 * lss_aniso_lambda[2],
|
||||
real_t(1.0) - Dplus0 * lss_aniso_lambda[0],
|
||||
real_t(1.0) - Dplus0 * lss_aniso_lambda[1],
|
||||
real_t(1.0) - Dplus0 * lss_aniso_lambda[2],
|
||||
};
|
||||
|
||||
//--------------------------------------------------------------------
|
||||
|
|
@ -360,7 +360,7 @@ int run( config_file& the_config )
|
|||
phi2.assign_function_of_grids_kdep([&](vec3_t<real_t> kvec, ccomplex_t pphi, ccomplex_t pphi2) {
|
||||
real_t k2 = kvec.norm_squared();
|
||||
real_t fac_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]);
|
||||
return pphi2 - (lss_aniso_sum_lambda * k2 + 4.0/3.0 * fac_aniso ) * pphi;
|
||||
return pphi2 - (lss_aniso_sum_lambda * k2 + real_t(4.0/3.0) * fac_aniso ) * pphi;
|
||||
}, phi, phi2);
|
||||
}
|
||||
|
||||
|
|
@ -651,7 +651,7 @@ int run( config_file& the_config )
|
|||
|
||||
tmp.FourierTransformBackward(false);
|
||||
tmp.assign_function_of_grids_r([&](auto ppsi, auto pgrad_psi, auto prho) {
|
||||
return vunit * std::real((std::conj(ppsi) * pgrad_psi - ppsi * std::conj(pgrad_psi)) / ccomplex_t(0.0, 2.0 / hbar)/(1.0+prho));
|
||||
return vunit * std::real((std::conj(ppsi) * pgrad_psi - ppsi * std::conj(pgrad_psi)) / ccomplex_t(0.0, 2.0 / hbar)/real_t(1.0+prho));
|
||||
}, psi, grad_psi, rho);
|
||||
|
||||
fluid_component fc = (idim==0)? fluid_component::vx : ((idim==1)? fluid_component::vy : fluid_component::vz );
|
||||
|
|
|
|||
|
|
@ -213,7 +213,7 @@ int main( int argc, char** argv )
|
|||
music::ilog << std::setw(32) << std::left << "OS/Kernel version" << " : " << kinfo.kernel << " version " << kinfo.major << "." << kinfo.minor << " build " << kinfo.build_number << std::endl;
|
||||
|
||||
// FFTW related infos
|
||||
music::ilog << std::setw(32) << std::left << "FFTW version" << " : " << fftw_version << std::endl;
|
||||
music::ilog << std::setw(32) << std::left << "FFTW version" << " : " << FFTW_API(version) << std::endl;
|
||||
music::ilog << std::setw(32) << std::left << "FFTW supports multi-threading" << " : " << (CONFIG::FFTW_threads_ok? "yes" : "no") << std::endl;
|
||||
music::ilog << std::setw(32) << std::left << "FFTW mode" << " : ";
|
||||
#if defined(FFTW_MODE_PATIENT)
|
||||
|
|
|
|||
|
|
@ -300,21 +300,15 @@ music_wnoise_generator<T>::music_wnoise_generator(/*const*/ music_wnoise_generat
|
|||
*rfine = new real_t[(size_t)rc.res_ * (size_t)rc.res_ * 2 * ((size_t)rc.res_ / 2 + 1)],
|
||||
*rcoarse = new real_t[(size_t)res_ * (size_t)res_ * 2 * ((size_t)res_ / 2 + 1)];
|
||||
|
||||
fftw_complex
|
||||
*ccoarse = reinterpret_cast<fftw_complex *>(rcoarse),
|
||||
*cfine = reinterpret_cast<fftw_complex *>(rfine);
|
||||
complex_t
|
||||
*ccoarse = reinterpret_cast<complex_t *>(rcoarse),
|
||||
*cfine = reinterpret_cast<complex_t *>(rfine);
|
||||
|
||||
int nx(rc.res_), ny(rc.res_), nz(rc.res_), nxc(res_), nyc(res_), nzc(res_);
|
||||
|
||||
#ifdef USE_SINGLEPRECISION
|
||||
fftwf_plan
|
||||
pf = fftwf_plan_dft_r2c_3d(nx, ny, nz, rfine, cfine, FFTW_ESTIMATE),
|
||||
ipc = fftwf_plan_dft_c2r_3d(nxc, nyc, nzc, ccoarse, rcoarse, FFTW_ESTIMATE);
|
||||
#else
|
||||
fftw_plan
|
||||
pf = fftw_plan_dft_r2c_3d(nx, ny, nz, rfine, cfine, FFTW_ESTIMATE),
|
||||
ipc = fftw_plan_dft_c2r_3d(nxc, nyc, nzc, ccoarse, rcoarse, FFTW_ESTIMATE);
|
||||
#endif
|
||||
FFTW_API(plan)
|
||||
pf = FFTW_API(plan_dft_r2c_3d)(nx, ny, nz, rfine, cfine, FFTW_ESTIMATE),
|
||||
ipc = FFTW_API(plan_dft_c2r_3d)(nxc, nyc, nzc, ccoarse, rcoarse, FFTW_ESTIMATE);
|
||||
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < nx; i++)
|
||||
|
|
@ -429,19 +423,11 @@ music_wnoise_generator<T>::music_wnoise_generator(music_wnoise_generator<T> &rc,
|
|||
nxc = lx[0] / 2, nyc = lx[1] / 2, nzc = lx[2] / 2;
|
||||
|
||||
real_t *rfine = new real_t[nx * ny * (nz + 2l)];
|
||||
fftw_complex *cfine = reinterpret_cast<fftw_complex *>(rfine);
|
||||
|
||||
|
||||
#ifdef USE_SINGLEPRECISION
|
||||
fftwf_plan
|
||||
pf = fftwf_plan_dft_r2c_3d(nx, ny, nz, rfine, cfine, FFTW_ESTIMATE),
|
||||
ipf = fftwf_plan_dft_c2r_3d(nx, ny, nz, cfine, rfine, FFTW_ESTIMATE);
|
||||
#else
|
||||
fftw_plan
|
||||
pf = fftw_plan_dft_r2c_3d(nx, ny, nz, rfine, cfine, FFTW_ESTIMATE),
|
||||
ipf = fftw_plan_dft_c2r_3d(nx, ny, nz, cfine, rfine, FFTW_ESTIMATE);
|
||||
#endif
|
||||
complex_t *cfine = reinterpret_cast<complex_t *>(rfine);
|
||||
|
||||
FFTW_API(plan)
|
||||
pf = FFTW_API(plan_dft_r2c_3d)(nx, ny, nz, rfine, cfine, FFTW_ESTIMATE),
|
||||
ipf = FFTW_API(plan_dft_c2r_3d)(nx, ny, nz, cfine, rfine, FFTW_ESTIMATE);
|
||||
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < (int)nx; i++)
|
||||
|
|
@ -454,14 +440,10 @@ music_wnoise_generator<T>::music_wnoise_generator(music_wnoise_generator<T> &rc,
|
|||
//this->free_all_mem(); // temporarily free memory, allocate again later
|
||||
|
||||
real_t *rcoarse = new real_t[nxc * nyc * (nzc + 2)];
|
||||
fftw_complex *ccoarse = reinterpret_cast<fftw_complex *>(rcoarse);
|
||||
complex_t *ccoarse = reinterpret_cast<complex_t *>(rcoarse);
|
||||
|
||||
|
||||
#ifdef USE_SINGLEPRECISION
|
||||
fftwf_plan pc = fftwf_plan_dft_r2c_3d(nxc, nyc, nzc, rcoarse, ccoarse, FFTW_ESTIMATE);
|
||||
#else
|
||||
fftw_plan pc = fftw_plan_dft_r2c_3d(nxc, nyc, nzc, rcoarse, ccoarse, FFTW_ESTIMATE);
|
||||
#endif
|
||||
FFTW_API(plan) pc = FFTW_API(plan_dft_r2c_3d)(nxc, nyc, nzc, rcoarse, ccoarse, FFTW_ESTIMATE);
|
||||
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < (int)nxc; i++)
|
||||
|
|
@ -824,5 +806,10 @@ void music_wnoise_generator<T>::print_allocated(void)
|
|||
music::ilog.Print(" -> %d of %d random number cubes currently allocated", ncount, ntot);
|
||||
}
|
||||
|
||||
#if defined(USE_PRECISION_FLOAT)
|
||||
template class music_wnoise_generator<float>;
|
||||
#elif defined(USE_PRECISION_DOUBLE)
|
||||
template class music_wnoise_generator<double>;
|
||||
#elif defined(USE_PRECISION_LONGDOUBLE)
|
||||
template class music_wnoise_generator<long double>;
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -377,15 +377,15 @@ public:
|
|||
auto argy = 0.5 * M_PI * kvec[1] / g.kny_[1];
|
||||
auto argz = 0.5 * M_PI * kvec[2] / g.kny_[2];
|
||||
|
||||
auto fx = sinc(argx);
|
||||
auto fx = real_t(sinc(argx));
|
||||
auto gx = ccomplex_t(0.0, dsinc(argx));
|
||||
auto fy = sinc(argy);
|
||||
auto fy = real_t(sinc(argy));
|
||||
auto gy = ccomplex_t(0.0, dsinc(argy));
|
||||
auto fz = sinc(argz);
|
||||
auto fz = real_t(sinc(argz));
|
||||
auto gz = ccomplex_t(0.0, dsinc(argz));
|
||||
|
||||
auto temp = (fx + sqrt3 * gx) * (fy + sqrt3 * gy) * (fz + sqrt3 * gz);
|
||||
auto magnitude = std::sqrt(1.0 - std::fabs(temp * temp));
|
||||
auto magnitude = real_t(std::sqrt(1.0 - std::fabs(temp * temp)));
|
||||
|
||||
auto y0(g0.kelem(i, j, k)), y1(g1.kelem(i, j, k)), y2(g2.kelem(i, j, k)), y3(g3.kelem(i, j, k)), y4(g4.kelem(i, j, k));
|
||||
|
||||
|
|
@ -519,16 +519,16 @@ public:
|
|||
auto argy = 0.5 * M_PI * kvec[1] / g.kny_[1];
|
||||
auto argz = 0.5 * M_PI * kvec[2] / g.kny_[2];
|
||||
|
||||
auto fx = sinc(argx);
|
||||
auto fx = real_t(sinc(argx));
|
||||
auto gx = ccomplex_t(0.0, dsinc(argx));
|
||||
auto fy = sinc(argy);
|
||||
auto fy = real_t(sinc(argy));
|
||||
auto gy = ccomplex_t(0.0, dsinc(argy));
|
||||
auto fz = sinc(argz);
|
||||
auto fz = real_t(sinc(argz));
|
||||
auto gz = ccomplex_t(0.0, dsinc(argz));
|
||||
|
||||
auto y1(g1.kelem(i, j, k)), y2(g2.kelem(i, j, k)), y3(g3.kelem(i, j, k)), y4(g4.kelem(i, j, k));
|
||||
|
||||
g0.kelem(i, j, k) += 3.0 * (y1 * gx * gy * fz + y2 * fx * gy * gz + y3 * gx * fy * gz) + sqrt27 * y4 * gx * gy * gz;
|
||||
g0.kelem(i, j, k) += real_t(3.0) * (y1 * gx * gy * fz + y2 * fx * gy * gz + y3 * gx * fy * gz) + sqrt27 * y4 * gx * gy * gz;
|
||||
|
||||
// fix the overall minus w.r.t. the monofonic noise definition
|
||||
g0.kelem(i, j, k ) *= -1.0;
|
||||
|
|
|
|||
|
|
@ -295,7 +295,7 @@ void output_convergence(
|
|||
for (std::size_t k = 0; k < phi_code.size(2); ++k) {
|
||||
std::size_t idx = phi_code.get_idx(i, j, k);
|
||||
auto kk = phi_code.get_k<real_t>(i, j, k);
|
||||
nabla_vini_mn.kelem(idx) = phi_code.kelem(idx) * (kk[m] * kk[n]);
|
||||
nabla_vini_mn.kelem(idx) = phi_code.kelem(idx) * real_t(kk[m] * kk[n]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -379,12 +379,12 @@ void output_convergence(
|
|||
for (std::size_t k = 0; k < phi.size(2); ++k) {
|
||||
auto kk = phi.get_k<real_t>(i, j, k);
|
||||
std::size_t idx = phi.get_idx(i, j, k);
|
||||
psi_1_tmp.kelem(idx) = ccomplex_t(0.0, 1.0) * (kk[idim] * phi.kelem(idx));
|
||||
psi_2_tmp.kelem(idx) = ccomplex_t(0.0, 1.0) * (kk[idim] * phi2.kelem(idx));
|
||||
psi_1_tmp.kelem(idx) = ccomplex_t(0.0, 1.0) * (real_t(kk[idim]) * phi.kelem(idx));
|
||||
psi_2_tmp.kelem(idx) = ccomplex_t(0.0, 1.0) * (real_t(kk[idim]) * phi2.kelem(idx));
|
||||
psi_3_tmp.kelem(idx) = ccomplex_t(0.0, 1.0) * (
|
||||
kk[idim] * phi3.kelem(idx) +
|
||||
kk[idimp] * A3[idimpp]->kelem(idx) -
|
||||
kk[idimpp] * A3[idimp]->kelem(idx)
|
||||
real_t(kk[idim]) * phi3.kelem(idx) +
|
||||
real_t(kk[idimp]) * A3[idimpp]->kelem(idx) -
|
||||
real_t(kk[idimpp]) * A3[idimp]->kelem(idx)
|
||||
);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in a new issue