mirror of
https://github.com/cosmo-sims/monofonIC.git
synced 2024-09-19 17:03:45 +02:00
added output to gadget2 particle files
This commit is contained in:
parent
9549f5f195
commit
17f2f7e8b6
6 changed files with 2141 additions and 377 deletions
11
example.conf
11
example.conf
|
@ -1,11 +1,18 @@
|
|||
[setup]
|
||||
LPTorder = 3
|
||||
GridRes = 128
|
||||
BoxLength = 300
|
||||
Dplus0 = 1.0
|
||||
BoxLength = 100
|
||||
zstart = 1.0
|
||||
H0 = 70.0
|
||||
Omega_m = 0.3
|
||||
Omega_L = 0.7
|
||||
Omega_b = 0.0455
|
||||
|
||||
[output]
|
||||
fname_hdf5 = output.hdf5
|
||||
fbase_analysis = output
|
||||
format = gadget2
|
||||
filename = ics_gadget.dat
|
||||
|
||||
[cosmology]
|
||||
transfer = eisenstein
|
||||
|
|
|
@ -1,4 +1,4 @@
|
|||
#pragma once
|
||||
#pragma once
|
||||
|
||||
#include <cmath>
|
||||
#include <array>
|
||||
|
@ -17,374 +17,395 @@ enum space_t
|
|||
template <typename data_t>
|
||||
class Grid_FFT
|
||||
{
|
||||
protected:
|
||||
protected:
|
||||
#if defined(USE_MPI)
|
||||
const MPI_Datatype MPI_data_t_type = (typeid(data_t) == typeid(double)) ? MPI_DOUBLE
|
||||
: (typeid(data_t) == typeid(float)) ? MPI_FLOAT
|
||||
: (typeid(data_t) == typeid(std::complex<float>)) ? MPI_COMPLEX
|
||||
: (typeid(data_t) == typeid(std::complex<double>)) ? MPI_DOUBLE_COMPLEX : MPI_INT;
|
||||
const MPI_Datatype MPI_data_t_type = (typeid(data_t) == typeid(double)) ? MPI_DOUBLE
|
||||
: (typeid(data_t) == typeid(float)) ? MPI_FLOAT
|
||||
: (typeid(data_t) == typeid(std::complex<float>)) ? MPI_COMPLEX
|
||||
: (typeid(data_t) == typeid(std::complex<double>)) ? MPI_DOUBLE_COMPLEX : MPI_INT;
|
||||
#endif
|
||||
public:
|
||||
std::array<size_t,3> n_, nhalf_;
|
||||
std::array<size_t,4> sizes_;
|
||||
size_t npr_, npc_;
|
||||
size_t ntot_;
|
||||
std::array<real_t,3> length_, kfac_, dx_;
|
||||
std::array<size_t, 3> n_, nhalf_;
|
||||
std::array<size_t, 4> sizes_;
|
||||
size_t npr_, npc_;
|
||||
size_t ntot_;
|
||||
std::array<real_t, 3> length_, kfac_, dx_;
|
||||
|
||||
space_t space_;
|
||||
data_t *data_;
|
||||
ccomplex_t *cdata_;
|
||||
space_t space_;
|
||||
data_t *data_;
|
||||
ccomplex_t *cdata_;
|
||||
|
||||
bounding_box<size_t> global_range_;
|
||||
bounding_box<size_t> global_range_;
|
||||
|
||||
fftw_plan_t plan_, iplan_;
|
||||
fftw_plan_t plan_, iplan_;
|
||||
|
||||
real_t fft_norm_fac_;
|
||||
real_t fft_norm_fac_;
|
||||
|
||||
ptrdiff_t local_0_start_, local_1_start_;
|
||||
ptrdiff_t local_0_size_, local_1_size_;
|
||||
ptrdiff_t local_0_start_, local_1_start_;
|
||||
ptrdiff_t local_0_size_, local_1_size_;
|
||||
|
||||
|
||||
|
||||
Grid_FFT(const std::array<size_t, 3> &N, const std::array<real_t, 3> &L, space_t initialspace = rspace_id)
|
||||
: n_(N), length_(L), space_(initialspace), data_(nullptr), cdata_(nullptr) //, RV_(*this), KV_(*this)
|
||||
{
|
||||
for(int i=0; i<3; ++i){
|
||||
kfac_[i] = 2.0 * M_PI / length_[i];
|
||||
dx_[i] = length_[i] / n_[i];
|
||||
}
|
||||
//invalidated = true;
|
||||
this->Setup();
|
||||
}
|
||||
|
||||
Grid_FFT(const Grid_FFT<data_t> &g)
|
||||
: n_(g.n_), length_(g.length_), space_(g.space_), data_(nullptr), cdata_(nullptr)
|
||||
{
|
||||
for (int i = 0; i < 3; ++i)
|
||||
{
|
||||
kfac_[i] = g.kfac_[i];
|
||||
dx_[i] = g.dx_[i];
|
||||
}
|
||||
//invalidated = true;
|
||||
this->Setup();
|
||||
|
||||
for (size_t i = 0; i < ntot_; ++i ){
|
||||
data_[i] = g.data_[i];
|
||||
}
|
||||
}
|
||||
|
||||
~Grid_FFT()
|
||||
{
|
||||
if( data_!=nullptr){
|
||||
fftw_free( data_ );
|
||||
}
|
||||
}
|
||||
|
||||
void Setup();
|
||||
|
||||
size_t size( size_t i ) const{ return sizes_[i]; }
|
||||
|
||||
void zero() {
|
||||
for( size_t i=0; i<ntot_; ++i ) data_[i] = 0.0;
|
||||
}
|
||||
|
||||
data_t &relem(size_t i, size_t j, size_t k)
|
||||
{
|
||||
size_t idx = (i * sizes_[1] + j) * sizes_[3] + k;
|
||||
return data_[idx];
|
||||
}
|
||||
|
||||
const data_t &relem(size_t i, size_t j, size_t k) const
|
||||
{
|
||||
size_t idx = (i * sizes_[1] + j) * sizes_[3] + k;
|
||||
return data_[idx];
|
||||
}
|
||||
|
||||
ccomplex_t &kelem(size_t i, size_t j, size_t k)
|
||||
{
|
||||
size_t idx = (i * sizes_[1] + j) * sizes_[3] + k;
|
||||
return cdata_[idx];
|
||||
}
|
||||
|
||||
const ccomplex_t &kelem(size_t i, size_t j, size_t k) const
|
||||
{
|
||||
size_t idx = (i * sizes_[1] + j) * sizes_[3] + k;
|
||||
return cdata_[idx];
|
||||
}
|
||||
|
||||
ccomplex_t &kelem(size_t idx) { return cdata_[idx]; }
|
||||
const ccomplex_t &kelem(size_t idx) const { return cdata_[idx]; }
|
||||
data_t &relem(size_t idx) { return data_[idx]; }
|
||||
const data_t &relem(size_t idx) const { return data_[idx]; }
|
||||
|
||||
size_t get_idx( size_t i, size_t j, size_t k ) const
|
||||
{
|
||||
return (i * sizes_[1] + j) * sizes_[3] + k;
|
||||
}
|
||||
|
||||
template <typename ft>
|
||||
vec3<ft> get_r(const size_t& i, const size_t& j, const size_t& k) const
|
||||
{
|
||||
vec3<ft> rr;
|
||||
|
||||
#if defined(USE_MPI)
|
||||
rr[0] = real_t(i + local_0_start_) * dx_[0];
|
||||
#else
|
||||
rr[0] = real_t(i) * dx_[0];
|
||||
#endif
|
||||
rr[1] = real_t(j) * dx_[1];
|
||||
rr[2] = real_t(k) * dx_[2];
|
||||
|
||||
return rr;
|
||||
}
|
||||
|
||||
template <typename ft>
|
||||
vec3<ft> get_k(const size_t &i, const size_t &j, const size_t &k) const
|
||||
{
|
||||
vec3<ft> kk;
|
||||
|
||||
#if defined(USE_MPI)
|
||||
auto ip = i+local_1_start_;
|
||||
kk[0] = (real_t(j) - real_t(j > nhalf_[0])*n_[0]) * kfac_[0];
|
||||
kk[1] = (real_t(ip) - real_t(ip > nhalf_[1])*n_[1]) * kfac_[1];
|
||||
#else
|
||||
kk[0] = (real_t(i) - real_t(i > nhalf_[0]) * n_[0]) * kfac_[0];
|
||||
kk[1] = (real_t(j) - real_t(j > nhalf_[1]) * n_[1]) * kfac_[1];
|
||||
#endif
|
||||
kk[2] = (real_t(k) - real_t(k>nhalf_[2])*n_[2])* kfac_[2];
|
||||
|
||||
return kk;
|
||||
}
|
||||
|
||||
template <typename functional>
|
||||
void apply_function_k(const functional &f)
|
||||
{
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
auto &elem = this->kelem(i, j, k);
|
||||
elem = f(elem);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <typename functional>
|
||||
void apply_function_r(const functional &f)
|
||||
{
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
auto &elem = this->relem(i, j, k);
|
||||
elem = f(elem);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double compute_2norm(void)
|
||||
{
|
||||
real_t sum1{0.0};
|
||||
#pragma omp parallel for reduction(+ : sum1)
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
const auto re = std::real(this->relem(i, j, k));
|
||||
const auto im = std::imag(this->relem(i, j, k));
|
||||
sum1 += re*re+im*im;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
sum1 /= sizes_[0] * sizes_[1] * sizes_[2];
|
||||
|
||||
return sum1;
|
||||
}
|
||||
|
||||
double std( void )
|
||||
{
|
||||
real_t sum1{0.0}, sum2{0.0};
|
||||
#pragma omp parallel for reduction(+:sum1,sum2)
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
const auto elem = std::real(this->relem(i, j, k));
|
||||
sum1 += elem;
|
||||
sum2 += elem*elem;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
sum1 /= sizes_[0] * sizes_[1] * sizes_[2];
|
||||
sum2 /= sizes_[0] * sizes_[1] * sizes_[2];
|
||||
|
||||
return std::sqrt( sum2 - sum1*sum1 );
|
||||
}
|
||||
double mean(void)
|
||||
{
|
||||
real_t sum1{0.0};
|
||||
#pragma omp parallel for reduction(+ : sum1)
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
const auto elem = std::real(this->relem(i, j, k));
|
||||
sum1 += elem;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
sum1 /= sizes_[0] * sizes_[1] * sizes_[2];
|
||||
|
||||
return sum1;
|
||||
}
|
||||
|
||||
template <typename functional, typename grid_t>
|
||||
void assign_function_of_grids_r(const functional &f, const grid_t &g)
|
||||
{
|
||||
assert(g.size(0) == size(0) && g.size(1) == size(1));// && g.size(2) == size(2) );
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
Grid_FFT(const std::array<size_t, 3> &N, const std::array<real_t, 3> &L, space_t initialspace = rspace_id)
|
||||
: n_(N), length_(L), space_(initialspace), data_(nullptr), cdata_(nullptr) //, RV_(*this), KV_(*this)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
for (int i = 0; i < 3; ++i)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
auto &elem = this->relem(i, j, k);
|
||||
const auto &elemg = g.relem(i, j, k);
|
||||
kfac_[i] = 2.0 * M_PI / length_[i];
|
||||
dx_[i] = length_[i] / n_[i];
|
||||
}
|
||||
//invalidated = true;
|
||||
this->Setup();
|
||||
}
|
||||
|
||||
elem = f(elemg);
|
||||
Grid_FFT(const Grid_FFT<data_t> &g)
|
||||
: n_(g.n_), length_(g.length_), space_(g.space_), data_(nullptr), cdata_(nullptr)
|
||||
{
|
||||
for (int i = 0; i < 3; ++i)
|
||||
{
|
||||
kfac_[i] = g.kfac_[i];
|
||||
dx_[i] = g.dx_[i];
|
||||
}
|
||||
//invalidated = true;
|
||||
this->Setup();
|
||||
|
||||
for (size_t i = 0; i < ntot_; ++i)
|
||||
{
|
||||
data_[i] = g.data_[i];
|
||||
}
|
||||
}
|
||||
|
||||
~Grid_FFT()
|
||||
{
|
||||
if (data_ != nullptr)
|
||||
{
|
||||
fftw_free(data_);
|
||||
}
|
||||
}
|
||||
|
||||
const Grid_FFT<data_t>* get_grid( size_t ilevel ) const { return this; }
|
||||
bool is_in_mask( size_t ilevel, size_t i, size_t j, size_t k ) const { return true; }
|
||||
bool is_refined( size_t ilevel, size_t i, size_t j, size_t k ) const { return false; }
|
||||
size_t levelmin() const {return 7;}
|
||||
size_t levelmax() const {return 7;}
|
||||
|
||||
void Setup();
|
||||
|
||||
size_t size(size_t i) const { return sizes_[i]; }
|
||||
|
||||
void zero()
|
||||
{
|
||||
for (size_t i = 0; i < ntot_; ++i)
|
||||
data_[i] = 0.0;
|
||||
}
|
||||
|
||||
data_t &relem(size_t i, size_t j, size_t k)
|
||||
{
|
||||
size_t idx = (i * sizes_[1] + j) * sizes_[3] + k;
|
||||
return data_[idx];
|
||||
}
|
||||
|
||||
const data_t &relem(size_t i, size_t j, size_t k) const
|
||||
{
|
||||
size_t idx = (i * sizes_[1] + j) * sizes_[3] + k;
|
||||
return data_[idx];
|
||||
}
|
||||
|
||||
ccomplex_t &kelem(size_t i, size_t j, size_t k)
|
||||
{
|
||||
size_t idx = (i * sizes_[1] + j) * sizes_[3] + k;
|
||||
return cdata_[idx];
|
||||
}
|
||||
|
||||
const ccomplex_t &kelem(size_t i, size_t j, size_t k) const
|
||||
{
|
||||
size_t idx = (i * sizes_[1] + j) * sizes_[3] + k;
|
||||
return cdata_[idx];
|
||||
}
|
||||
|
||||
ccomplex_t &kelem(size_t idx) { return cdata_[idx]; }
|
||||
const ccomplex_t &kelem(size_t idx) const { return cdata_[idx]; }
|
||||
data_t &relem(size_t idx) { return data_[idx]; }
|
||||
const data_t &relem(size_t idx) const { return data_[idx]; }
|
||||
|
||||
size_t get_idx(size_t i, size_t j, size_t k) const
|
||||
{
|
||||
return (i * sizes_[1] + j) * sizes_[3] + k;
|
||||
}
|
||||
|
||||
template <typename ft>
|
||||
vec3<ft> get_r(const size_t &i, const size_t &j, const size_t &k) const
|
||||
{
|
||||
vec3<ft> rr;
|
||||
|
||||
#if defined(USE_MPI)
|
||||
rr[0] = real_t(i + local_0_start_) * dx_[0];
|
||||
#else
|
||||
rr[0] = real_t(i) * dx_[0];
|
||||
#endif
|
||||
rr[1] = real_t(j) * dx_[1];
|
||||
rr[2] = real_t(k) * dx_[2];
|
||||
|
||||
return rr;
|
||||
}
|
||||
|
||||
void cell_pos( int ilevel, size_t i, size_t j, size_t k, double* x ) const {
|
||||
|
||||
x[0] = double(i)/size(0);
|
||||
x[1] = double(j)/size(1);
|
||||
x[2] = double(k)/size(2);
|
||||
}
|
||||
|
||||
size_t count_leaf_cells( int, int ) const {
|
||||
return n_[0]*n_[1]*n_[2];
|
||||
}
|
||||
|
||||
template <typename ft>
|
||||
vec3<ft> get_k(const size_t &i, const size_t &j, const size_t &k) const
|
||||
{
|
||||
vec3<ft> kk;
|
||||
|
||||
#if defined(USE_MPI)
|
||||
auto ip = i + local_1_start_;
|
||||
kk[0] = (real_t(j) - real_t(j > nhalf_[0]) * n_[0]) * kfac_[0];
|
||||
kk[1] = (real_t(ip) - real_t(ip > nhalf_[1]) * n_[1]) * kfac_[1];
|
||||
#else
|
||||
kk[0] = (real_t(i) - real_t(i > nhalf_[0]) * n_[0]) * kfac_[0];
|
||||
kk[1] = (real_t(j) - real_t(j > nhalf_[1]) * n_[1]) * kfac_[1];
|
||||
#endif
|
||||
kk[2] = (real_t(k) - real_t(k > nhalf_[2]) * n_[2]) * kfac_[2];
|
||||
|
||||
return kk;
|
||||
}
|
||||
|
||||
template <typename functional>
|
||||
void apply_function_k(const functional &f)
|
||||
{
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
auto &elem = this->kelem(i, j, k);
|
||||
elem = f(elem);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template <typename functional>
|
||||
void apply_function_r(const functional &f)
|
||||
{
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
auto &elem = this->relem(i, j, k);
|
||||
elem = f(elem);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <typename functional, typename grid1_t, typename grid2_t>
|
||||
void assign_function_of_grids_r(const functional &f, const grid1_t &g1, const grid2_t &g2)
|
||||
{
|
||||
assert(g1.size(0) == size(0) && g1.size(1) == size(1));// && g1.size(2) == size(2));
|
||||
assert(g2.size(0) == size(0) && g2.size(1) == size(1));// && g2.size(2) == size(2));
|
||||
double compute_2norm(void)
|
||||
{
|
||||
real_t sum1{0.0};
|
||||
#pragma omp parallel for reduction(+ \
|
||||
: sum1)
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
const auto re = std::real(this->relem(i, j, k));
|
||||
const auto im = std::imag(this->relem(i, j, k));
|
||||
sum1 += re * re + im * im;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
sum1 /= sizes_[0] * sizes_[1] * sizes_[2];
|
||||
|
||||
return sum1;
|
||||
}
|
||||
|
||||
double std(void)
|
||||
{
|
||||
real_t sum1{0.0}, sum2{0.0};
|
||||
#pragma omp parallel for reduction(+ \
|
||||
: sum1, sum2)
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
const auto elem = std::real(this->relem(i, j, k));
|
||||
sum1 += elem;
|
||||
sum2 += elem * elem;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
sum1 /= sizes_[0] * sizes_[1] * sizes_[2];
|
||||
sum2 /= sizes_[0] * sizes_[1] * sizes_[2];
|
||||
|
||||
return std::sqrt(sum2 - sum1 * sum1);
|
||||
}
|
||||
double mean(void)
|
||||
{
|
||||
real_t sum1{0.0};
|
||||
#pragma omp parallel for reduction(+ \
|
||||
: sum1)
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
const auto elem = std::real(this->relem(i, j, k));
|
||||
sum1 += elem;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
sum1 /= sizes_[0] * sizes_[1] * sizes_[2];
|
||||
|
||||
return sum1;
|
||||
}
|
||||
|
||||
template <typename functional, typename grid_t>
|
||||
void assign_function_of_grids_r(const functional &f, const grid_t &g)
|
||||
{
|
||||
assert(g.size(0) == size(0) && g.size(1) == size(1)); // && g.size(2) == size(2) );
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
//auto idx = this->get_idx(i,j,k);
|
||||
auto &elem = this->relem(i,j,k);
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
auto &elem = this->relem(i, j, k);
|
||||
const auto &elemg = g.relem(i, j, k);
|
||||
|
||||
const auto &elemg1 = g1.relem(i,j,k);
|
||||
const auto &elemg2 = g2.relem(i,j,k);
|
||||
elem = f(elemg);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
elem = f(elemg1,elemg2);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <typename functional, typename grid1_t, typename grid2_t, typename grid3_t>
|
||||
void assign_function_of_grids_r(const functional &f, const grid1_t &g1, const grid2_t &g2, const grid3_t &g3)
|
||||
{
|
||||
assert(g1.size(0) == size(0) && g1.size(1) == size(1));// && g1.size(2) == size(2));
|
||||
assert(g2.size(0) == size(0) && g2.size(1) == size(1));// && g2.size(2) == size(2));
|
||||
assert(g3.size(0) == size(0) && g3.size(1) == size(1));// && g3.size(2) == size(2));
|
||||
template <typename functional, typename grid1_t, typename grid2_t>
|
||||
void assign_function_of_grids_r(const functional &f, const grid1_t &g1, const grid2_t &g2)
|
||||
{
|
||||
assert(g1.size(0) == size(0) && g1.size(1) == size(1)); // && g1.size(2) == size(2));
|
||||
assert(g2.size(0) == size(0) && g2.size(1) == size(1)); // && g2.size(2) == size(2));
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
//auto idx = this->get_idx(i,j,k);
|
||||
auto &elem = this->relem(i, j, k);
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
//auto idx = this->get_idx(i,j,k);
|
||||
auto &elem = this->relem(i, j, k);
|
||||
|
||||
const auto &elemg1 = g1.relem(i, j, k);
|
||||
const auto &elemg2 = g2.relem(i, j, k);
|
||||
const auto &elemg3 = g3.relem(i, j, k);
|
||||
const auto &elemg1 = g1.relem(i, j, k);
|
||||
const auto &elemg2 = g2.relem(i, j, k);
|
||||
|
||||
elem = f(elemg1, elemg2, elemg3);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
elem = f(elemg1, elemg2);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template< typename functional >
|
||||
void apply_function_k_dep( const functional& f )
|
||||
{
|
||||
#pragma omp parallel for
|
||||
for( size_t i = 0; i < sizes_[0]; ++i )
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
auto& elem = this->kelem(i,j,k);
|
||||
elem = f(elem,this->get_k<real_t>(i,j,k));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
template <typename functional, typename grid1_t, typename grid2_t, typename grid3_t>
|
||||
void assign_function_of_grids_r(const functional &f, const grid1_t &g1, const grid2_t &g2, const grid3_t &g3)
|
||||
{
|
||||
assert(g1.size(0) == size(0) && g1.size(1) == size(1)); // && g1.size(2) == size(2));
|
||||
assert(g2.size(0) == size(0) && g2.size(1) == size(1)); // && g2.size(2) == size(2));
|
||||
assert(g3.size(0) == size(0) && g3.size(1) == size(1)); // && g3.size(2) == size(2));
|
||||
|
||||
template <typename functional>
|
||||
void apply_function_r_dep( const functional& f)
|
||||
{
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
auto &elem = this->relem(i, j, k);
|
||||
elem = f(elem, this->get_r<real_t>(i, j, k));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
//auto idx = this->get_idx(i,j,k);
|
||||
auto &elem = this->relem(i, j, k);
|
||||
|
||||
void FourierTransformBackward(bool do_transform = true);
|
||||
const auto &elemg1 = g1.relem(i, j, k);
|
||||
const auto &elemg2 = g2.relem(i, j, k);
|
||||
const auto &elemg3 = g3.relem(i, j, k);
|
||||
|
||||
void FourierTransformForward(bool do_transform = true);
|
||||
elem = f(elemg1, elemg2, elemg3);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void ApplyNorm(void);
|
||||
template <typename functional>
|
||||
void apply_function_k_dep(const functional &f)
|
||||
{
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
auto &elem = this->kelem(i, j, k);
|
||||
elem = f(elem, this->get_k<real_t>(i, j, k));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void FillRandomReal( unsigned long int seed = 123456ul );
|
||||
template <typename functional>
|
||||
void apply_function_r_dep(const functional &f)
|
||||
{
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
auto &elem = this->relem(i, j, k);
|
||||
elem = f(elem, this->get_r<real_t>(i, j, k));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void Write_to_HDF5(std::string fname, std::string datasetname);
|
||||
void FourierTransformBackward(bool do_transform = true);
|
||||
|
||||
void ComputePowerSpectrum(std::vector<double> &bin_k, std::vector<double> &bin_P, std::vector<double> &bin_eP, int nbins);
|
||||
void FourierTransformForward(bool do_transform = true);
|
||||
|
||||
void Compute_PDF(std::string ofname, int nbins = 1000, double scale = 1.0, double rhomin = 1e-3, double rhomax = 1e3);
|
||||
void ApplyNorm(void);
|
||||
|
||||
void zero_DC_mode(void)
|
||||
{
|
||||
void FillRandomReal(unsigned long int seed = 123456ul);
|
||||
|
||||
void Write_to_HDF5(std::string fname, std::string datasetname);
|
||||
|
||||
void ComputePowerSpectrum(std::vector<double> &bin_k, std::vector<double> &bin_P, std::vector<double> &bin_eP, int nbins);
|
||||
|
||||
void Compute_PDF(std::string ofname, int nbins = 1000, double scale = 1.0, double rhomin = 1e-3, double rhomax = 1e3);
|
||||
|
||||
void zero_DC_mode(void)
|
||||
{
|
||||
#ifdef USE_MPI
|
||||
if (CONFIG::MPI_task_rank == 0)
|
||||
if (CONFIG::MPI_task_rank == 0)
|
||||
#endif
|
||||
cdata_[0] = (data_t)0.0;
|
||||
}
|
||||
cdata_[0] = (data_t)0.0;
|
||||
}
|
||||
};
|
||||
|
||||
template <typename data_t>
|
||||
|
|
169
include/output_plugin.hh
Normal file
169
include/output_plugin.hh
Normal file
|
@ -0,0 +1,169 @@
|
|||
/*
|
||||
|
||||
output.hh - This file is part of MUSIC -
|
||||
a code to generate multi-scale initial conditions
|
||||
for cosmological simulations
|
||||
|
||||
Copyright (C) 2010 Oliver Hahn
|
||||
|
||||
*/
|
||||
|
||||
#ifndef __OUTPUT_HH
|
||||
#define __OUTPUT_HH
|
||||
|
||||
#include <string>
|
||||
#include <map>
|
||||
|
||||
#include "general.hh"
|
||||
#include "grid_fft.hh"
|
||||
#include "config_file.hh"
|
||||
|
||||
/*!
|
||||
* @class output_plugin
|
||||
* @brief abstract base class for output plug-ins
|
||||
*
|
||||
* This class provides the abstract base class for all output plug-ins.
|
||||
* All output plug-ins need to derive from it and implement the purely
|
||||
* virtual member functions.
|
||||
*/
|
||||
class output_plugin
|
||||
{
|
||||
|
||||
public:
|
||||
using grid_hierarchy = Grid_FFT<real_t>;
|
||||
|
||||
protected:
|
||||
|
||||
//! reference to the ConfigFile object that holds all configuration options
|
||||
ConfigFile& cf_;
|
||||
|
||||
//! output file or directory name
|
||||
std::string fname_;
|
||||
|
||||
//! minimum refinement level
|
||||
// unsigned levelmin_;
|
||||
|
||||
//! maximum refinement level
|
||||
// unsigned levelmax_;
|
||||
|
||||
std::vector<unsigned>
|
||||
offx_, //!< vector describing the x-offset of each level
|
||||
offy_, //!< vector describing the y-offset of each level
|
||||
offz_, //!< vector describing the z-offset of each level
|
||||
sizex_, //!< vector describing the extent in x of each level
|
||||
sizey_, //!< vector describing the extent in y of each level
|
||||
sizez_; //!< vector describing the extent in z of each level
|
||||
|
||||
//! quick access function to query properties of the refinement grid from the configuration options
|
||||
/*! @param name name of the config property
|
||||
* @param icomp component index (0=x, 1=y, 2=z)
|
||||
* @param oit output iterator (e.g. std::back_inserter for vectors)
|
||||
*/
|
||||
template< typename output_iterator >
|
||||
void query_grid_prop( std::string name, int icomp, output_iterator oit )
|
||||
{
|
||||
char str[128];
|
||||
//for( unsigned i=levelmin_; i<=levelmax_; ++i )
|
||||
unsigned i=0;
|
||||
{
|
||||
sprintf( str, "%s(%u,%d)", name.c_str(), i, icomp );
|
||||
*oit = 0; //cf_.GetValue<unsigned>( "setup", str );
|
||||
++oit;
|
||||
}
|
||||
}
|
||||
|
||||
public:
|
||||
|
||||
//! constructor
|
||||
explicit output_plugin( ConfigFile& cf )
|
||||
: cf_(cf)
|
||||
{
|
||||
fname_ = cf.GetValue<std::string>("output","filename");
|
||||
// levelmin_ = cf.GetValue<unsigned>( "setup", "levelmin" );
|
||||
// levelmax_ = cf.GetValue<unsigned>( "setup", "levelmax" );
|
||||
|
||||
query_grid_prop( "offset", 0, std::back_inserter(offx_) );
|
||||
query_grid_prop( "offset", 1, std::back_inserter(offy_) );
|
||||
query_grid_prop( "offset", 2, std::back_inserter(offz_) );
|
||||
|
||||
query_grid_prop( "size", 0, std::back_inserter(sizex_) );
|
||||
query_grid_prop( "size", 1, std::back_inserter(sizey_) );
|
||||
query_grid_prop( "size", 2, std::back_inserter(sizez_) );
|
||||
}
|
||||
|
||||
//! destructor
|
||||
virtual ~output_plugin()
|
||||
{ }
|
||||
|
||||
//! purely virtual prototype to write the masses for each dark matter particle
|
||||
virtual void write_dm_mass( const grid_hierarchy& gh ) = 0;
|
||||
|
||||
//! purely virtual prototype to write the dark matter density field
|
||||
virtual void write_dm_density( const grid_hierarchy& gh ) = 0;
|
||||
|
||||
//! purely virtual prototype to write the dark matter gravitational potential (from which displacements are computed in 1LPT)
|
||||
virtual void write_dm_potential( const grid_hierarchy& gh ) = 0;
|
||||
|
||||
//! purely virtual prototype to write dark matter particle velocities
|
||||
virtual void write_dm_velocity( int coord, const grid_hierarchy& gh ) = 0;
|
||||
|
||||
//! purely virtual prototype to write dark matter particle positions
|
||||
virtual void write_dm_position( int coord, const grid_hierarchy& gh ) = 0;
|
||||
|
||||
//! purely virtual prototype to write the baryon velocities
|
||||
virtual void write_gas_velocity( int coord, const grid_hierarchy& gh ) = 0;
|
||||
|
||||
//! purely virtual prototype to write the baryon coordinates
|
||||
virtual void write_gas_position( int coord, const grid_hierarchy& gh ) = 0;
|
||||
|
||||
//! purely virtual prototype to write the baryon density field
|
||||
virtual void write_gas_density( const grid_hierarchy& gh ) = 0;
|
||||
|
||||
//! purely virtual prototype to write the baryon gravitational potential (from which displacements are computed in 1LPT)
|
||||
virtual void write_gas_potential( const grid_hierarchy& gh ) = 0;
|
||||
|
||||
//! purely virtual prototype for all things to be done at the very end
|
||||
virtual void finalize( void ) = 0;
|
||||
};
|
||||
|
||||
/*!
|
||||
* @brief implements abstract factory design pattern for output plug-ins
|
||||
*/
|
||||
struct output_plugin_creator
|
||||
{
|
||||
//! create an instance of a plug-in
|
||||
virtual output_plugin * create( ConfigFile& cf ) const = 0;
|
||||
|
||||
//! destroy an instance of a plug-in
|
||||
virtual ~output_plugin_creator() { }
|
||||
};
|
||||
|
||||
//! maps the name of a plug-in to a pointer of the factory pattern
|
||||
std::map< std::string, output_plugin_creator *>& get_output_plugin_map();
|
||||
|
||||
//! print a list of all registered output plug-ins
|
||||
void print_output_plugins();
|
||||
|
||||
/*!
|
||||
* @brief concrete factory pattern for output plug-ins
|
||||
*/
|
||||
template< class Derived >
|
||||
struct output_plugin_creator_concrete : public output_plugin_creator
|
||||
{
|
||||
//! register the plug-in by its name
|
||||
output_plugin_creator_concrete( const std::string& plugin_name )
|
||||
{
|
||||
get_output_plugin_map()[ plugin_name ] = this;
|
||||
}
|
||||
|
||||
//! create an instance of the plug-in
|
||||
output_plugin * create( ConfigFile& cf ) const
|
||||
{
|
||||
return new Derived( cf );
|
||||
}
|
||||
};
|
||||
|
||||
//! failsafe version to select the output plug-in
|
||||
output_plugin *select_output_plugin( ConfigFile& cf );
|
||||
|
||||
#endif // __OUTPUT_HH
|
186
src/main.cc
186
src/main.cc
|
@ -11,6 +11,7 @@
|
|||
|
||||
#include <transfer_function_plugin.hh>
|
||||
#include <random_plugin.hh>
|
||||
#include <output_plugin.hh>
|
||||
#include <cosmology_calculator.hh>
|
||||
|
||||
namespace CONFIG{
|
||||
|
@ -24,6 +25,7 @@ bool FFTW_threads_ok = false;
|
|||
|
||||
RNG_plugin *the_random_number_generator;
|
||||
TransferFunction_plugin *the_transfer_function;
|
||||
output_plugin *the_output_plugin;
|
||||
|
||||
int main( int argc, char** argv )
|
||||
{
|
||||
|
@ -66,7 +68,7 @@ int main( int argc, char** argv )
|
|||
// print_region_generator_plugins();
|
||||
print_TransferFunction_plugins();
|
||||
// print_RNG_plugins();
|
||||
// print_output_plugins();
|
||||
print_output_plugins();
|
||||
|
||||
csoca::elog << "In order to run, you need to specify a parameter file!" << std::endl;
|
||||
exit(0);
|
||||
|
@ -79,40 +81,41 @@ int main( int argc, char** argv )
|
|||
|
||||
const size_t ngrid = the_config.GetValue<size_t>("setup", "GridRes");
|
||||
const real_t boxlen = the_config.GetValue<double>("setup", "BoxLength");
|
||||
const real_t zstart = the_config.GetValue<double>("setup", "zstart");
|
||||
const int LPTorder = the_config.GetValueSafe<double>("setup","LPTorder",100);
|
||||
const real_t astart = 1.0/(1.0+zstart);
|
||||
const real_t volfac(std::pow(boxlen / ngrid / 2.0 / M_PI, 1.5));
|
||||
const real_t phifac = 1.0 / boxlen / boxlen; // to have potential in box units
|
||||
const real_t deriv_fac = 1.0 ;//boxlen;
|
||||
|
||||
real_t Dplus0 = the_config.GetValue<real_t>("setup", "Dplus0");
|
||||
// real_t Dplus0 = the_config.GetValue<real_t>("setup", "Dplus0");
|
||||
// real_t Ddot0 = 1.0;
|
||||
|
||||
const bool bDoFixing = false;
|
||||
|
||||
|
||||
|
||||
//...
|
||||
const std::string fname_hdf5 = the_config.GetValueSafe<std::string>("output", "fname_hdf5", "output.hdf5");
|
||||
//////////////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
std::unique_ptr<CosmologyCalculator>
|
||||
the_cosmo_calc;
|
||||
|
||||
|
||||
try
|
||||
{
|
||||
the_random_number_generator = select_RNG_plugin(the_config);
|
||||
the_transfer_function = select_TransferFunction_plugin(the_config);
|
||||
the_output_plugin = select_output_plugin(the_config);
|
||||
|
||||
the_cosmo_calc = std::make_unique<CosmologyCalculator>(the_config, the_transfer_function);
|
||||
|
||||
//double pnorm = the_cosmo_calc->ComputePNorm();
|
||||
//Dplus = the_cosmo_calc->CalcGrowthFactor(astart) / the_cosmo_calc->CalcGrowthFactor(1.0);
|
||||
|
||||
csoca::ilog << "power spectrum is output for D+ =" << Dplus0 << std::endl;
|
||||
// double pnorm = the_cosmo_calc->ComputePNorm();
|
||||
// csoca::ilog << "power spectrum is output for D+ =" << Dplus0 << std::endl;
|
||||
//csoca::ilog << "power spectrum normalisation is " << pnorm << std::endl;
|
||||
//csoca::ilog << "power spectrum normalisation is " << pnorm*Dplus*Dplus << std::endl;
|
||||
|
||||
// write power spectrum to a file
|
||||
std::ofstream ofs("input_powerspec.txt");
|
||||
for( double k=1e-4; k<1e4; k*=1.1 ){
|
||||
ofs << std::setw(16) << k
|
||||
<< std::setw(16) << std::pow(the_cosmo_calc->GetAmplitude(k, total) * Dplus0, 2.0)
|
||||
<< std::setw(16) << std::pow(the_cosmo_calc->GetAmplitude(k, total), 2.0)
|
||||
<< std::endl;
|
||||
}
|
||||
|
||||
}catch(...){
|
||||
csoca::elog << "Problem during initialisation. See error(s) above. Exiting..." << std::endl;
|
||||
#if defined(USE_MPI)
|
||||
|
@ -120,6 +123,30 @@ int main( int argc, char** argv )
|
|||
#endif
|
||||
return 1;
|
||||
}
|
||||
const real_t Dplus0 = the_cosmo_calc->CalcGrowthFactor(astart) / the_cosmo_calc->CalcGrowthFactor(1.0);
|
||||
const real_t vfac = the_cosmo_calc->CalcVFact(astart);
|
||||
|
||||
{
|
||||
// write power spectrum to a file
|
||||
std::ofstream ofs("input_powerspec.txt");
|
||||
for( double k=1e-4; k<1e4; k*=1.1 ){
|
||||
ofs << std::setw(16) << k
|
||||
<< std::setw(16) << std::pow(the_cosmo_calc->GetAmplitude(k, total) * Dplus0, 2.0)
|
||||
<< std::setw(16) << std::pow(the_cosmo_calc->GetAmplitude(k, total), 2.0)
|
||||
<< std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
// compute growth factors of the respective orders
|
||||
const double g1 = -Dplus0;
|
||||
const double g2 = -3.0/7.0*Dplus0*Dplus0;
|
||||
const double g3a = -1.0/3.0*Dplus0*Dplus0*Dplus0;
|
||||
const double g3b = 10.0/21.*Dplus0*Dplus0*Dplus0;
|
||||
|
||||
const double vfac1 = vfac;
|
||||
const double vfac2 = 2*vfac1;
|
||||
const double vfac3 = 3*vfac1;
|
||||
|
||||
//--------------------------------------------------------------------
|
||||
// Create arrays
|
||||
Grid_FFT<real_t> phi({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
|
@ -136,6 +163,8 @@ int main( int argc, char** argv )
|
|||
|
||||
phi.apply_function_k_dep([&](auto x, auto k) -> ccomplex_t {
|
||||
real_t kmod = k.norm();
|
||||
if( bDoFixing ) x = x / std::abs(x); //std::exp(ccomplex_t(0, iphase * PhaseRotation));
|
||||
else x = x;
|
||||
ccomplex_t delta = x * the_cosmo_calc->GetAmplitude(kmod, total);
|
||||
return -delta / (kmod * kmod) * phifac / volfac;
|
||||
});
|
||||
|
@ -295,16 +324,29 @@ int main( int argc, char** argv )
|
|||
phi3b.zero_DC_mode();
|
||||
|
||||
///////////////////////////////////////////////////////////////////////
|
||||
|
||||
// we store the densities here if we compute them
|
||||
Grid_FFT<real_t> &delta = phi_xx;
|
||||
Grid_FFT<real_t> &delta2 = phi_xy;
|
||||
Grid_FFT<real_t> &delta3a = phi_xz;
|
||||
Grid_FFT<real_t> &delta3b = phi_yy;
|
||||
Grid_FFT<real_t> &delta3 = phi_yz;
|
||||
|
||||
delta.FourierTransformForward(false);
|
||||
delta2.FourierTransformForward(false);
|
||||
delta3a.FourierTransformForward(false);
|
||||
delta3b.FourierTransformForward(false);
|
||||
// we store displacements and velocities here if we compute them
|
||||
Grid_FFT<real_t> &Psix = phi_xx;
|
||||
Grid_FFT<real_t> &Psiy = phi_xy;
|
||||
Grid_FFT<real_t> &Psiz = phi_xz;
|
||||
Grid_FFT<real_t> &Vx = phi_yy;
|
||||
Grid_FFT<real_t> &Vy = phi_yz;
|
||||
Grid_FFT<real_t> &Vz = phi_zz;
|
||||
|
||||
const bool compute_densities = false;
|
||||
|
||||
phi_xx.FourierTransformForward(false);
|
||||
phi_xy.FourierTransformForward(false);
|
||||
phi_xz.FourierTransformForward(false);
|
||||
phi_yy.FourierTransformForward(false);
|
||||
phi_yz.FourierTransformForward(false);
|
||||
phi_zz.FourierTransformForward(false);
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < phi.size(0); ++i)
|
||||
|
@ -317,36 +359,94 @@ int main( int argc, char** argv )
|
|||
size_t idx = phi.get_idx(i,j,k);
|
||||
auto laplace = -kk.norm_squared();
|
||||
|
||||
delta.kelem(idx) = laplace * phi.kelem(idx) / phifac;
|
||||
delta2.kelem(idx) = laplace * phi2.kelem(idx) / phifac;
|
||||
delta3a.kelem(idx) = laplace * phi3a.kelem(idx) / phifac;
|
||||
delta3b.kelem(idx) = laplace * phi3b.kelem(idx) / phifac;
|
||||
// scale potentials with respective order growth factors
|
||||
phi.kelem(idx) *= g1;
|
||||
phi2.kelem(idx) *= g2;
|
||||
phi3a.kelem(idx) *= g3a;
|
||||
phi3b.kelem(idx) *= g3b;
|
||||
|
||||
if( compute_densities ){
|
||||
// compute densities associated to respective potentials as well
|
||||
delta.kelem(idx) = laplace * phi.kelem(idx) / phifac;
|
||||
delta2.kelem(idx) = laplace * phi2.kelem(idx) / phifac;
|
||||
delta3a.kelem(idx) = laplace * phi3a.kelem(idx) / phifac;
|
||||
delta3b.kelem(idx) = laplace * phi3b.kelem(idx) / phifac;
|
||||
delta3.kelem(idx) = delta3a.kelem(idx) + delta3b.kelem(idx);
|
||||
}else{
|
||||
auto phitot = phi.kelem(idx) + ((LPTorder>1)?phi2.kelem(idx):0.0) + ((LPTorder>2)? phi3a.kelem(idx) + phi3b.kelem(idx) : 0.0);
|
||||
auto phitot_v = vfac1 * phi.kelem(idx) + ((LPTorder>1)? vfac2 * phi2.kelem(idx) : 0.0) + ((LPTorder>2)? vfac3 * (phi3a.kelem(idx) + phi3b.kelem(idx)) : 0.0);
|
||||
|
||||
Psix.kelem(idx) = ccomplex_t(0.0,1.0) * kk[0]* boxlen * ( phitot );
|
||||
Psiy.kelem(idx) = ccomplex_t(0.0,1.0) * kk[1]* boxlen * ( phitot );
|
||||
Psiz.kelem(idx) = ccomplex_t(0.0,1.0) * kk[2]* boxlen * ( phitot );
|
||||
|
||||
Vx.kelem(idx) = ccomplex_t(0.0,1.0) * kk[0]* boxlen * ( phitot_v );
|
||||
Vy.kelem(idx) = ccomplex_t(0.0,1.0) * kk[1]* boxlen * ( phitot_v );
|
||||
Vz.kelem(idx) = ccomplex_t(0.0,1.0) * kk[2]* boxlen * ( phitot_v );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
///////////////////////////////////////////////////////////////////////
|
||||
phi.FourierTransformBackward();
|
||||
phi2.FourierTransformBackward();
|
||||
phi3a.FourierTransformBackward();
|
||||
phi3b.FourierTransformBackward();
|
||||
|
||||
delta.FourierTransformBackward();
|
||||
delta2.FourierTransformBackward();
|
||||
delta3a.FourierTransformBackward();
|
||||
delta3b.FourierTransformBackward();
|
||||
|
||||
|
||||
|
||||
//... write output .....
|
||||
unlink(fname_hdf5.c_str());
|
||||
phi.Write_to_HDF5(fname_hdf5, "phi");
|
||||
phi2.Write_to_HDF5(fname_hdf5, "phi2");
|
||||
phi3a.Write_to_HDF5(fname_hdf5, "phi3a");
|
||||
phi3b.Write_to_HDF5(fname_hdf5, "phi3b");
|
||||
delta.Write_to_HDF5(fname_hdf5, "delta");
|
||||
delta2.Write_to_HDF5(fname_hdf5, "delta2");
|
||||
delta3a.Write_to_HDF5(fname_hdf5, "delta3a");
|
||||
delta3b.Write_to_HDF5(fname_hdf5, "delta3b");
|
||||
|
||||
|
||||
if( compute_densities ){
|
||||
phi.FourierTransformBackward();
|
||||
phi2.FourierTransformBackward();
|
||||
phi3a.FourierTransformBackward();
|
||||
phi3b.FourierTransformBackward();
|
||||
|
||||
delta.FourierTransformBackward();
|
||||
delta2.FourierTransformBackward();
|
||||
delta3a.FourierTransformBackward();
|
||||
delta3b.FourierTransformBackward();
|
||||
delta3.FourierTransformBackward();
|
||||
|
||||
unlink(fname_hdf5.c_str());
|
||||
phi.Write_to_HDF5(fname_hdf5, "phi");
|
||||
phi2.Write_to_HDF5(fname_hdf5, "phi2");
|
||||
phi3a.Write_to_HDF5(fname_hdf5, "phi3a");
|
||||
phi3b.Write_to_HDF5(fname_hdf5, "phi3b");
|
||||
|
||||
delta.Write_to_HDF5(fname_hdf5, "delta");
|
||||
delta2.Write_to_HDF5(fname_hdf5, "delta2");
|
||||
delta3a.Write_to_HDF5(fname_hdf5, "delta3a");
|
||||
delta3b.Write_to_HDF5(fname_hdf5, "delta3b");
|
||||
delta3.Write_to_HDF5(fname_hdf5, "delta3");
|
||||
}else{
|
||||
Psix.FourierTransformBackward();
|
||||
Psiy.FourierTransformBackward();
|
||||
Psiz.FourierTransformBackward();
|
||||
Vx.FourierTransformBackward();
|
||||
Vy.FourierTransformBackward();
|
||||
Vz.FourierTransformBackward();
|
||||
|
||||
// Psix.Write_to_HDF5(fname_hdf5, "Psix");
|
||||
// Psiy.Write_to_HDF5(fname_hdf5, "Psiy");
|
||||
// Psiz.Write_to_HDF5(fname_hdf5, "Psiz");
|
||||
// Vx.Write_to_HDF5(fname_hdf5, "Vx");
|
||||
// Vy.Write_to_HDF5(fname_hdf5, "Vy");
|
||||
// Vz.Write_to_HDF5(fname_hdf5, "Vz");
|
||||
|
||||
the_output_plugin->write_dm_mass(Psix);
|
||||
the_output_plugin->write_dm_density(Psix);
|
||||
|
||||
the_output_plugin->write_dm_position(0, Psix );
|
||||
the_output_plugin->write_dm_position(1, Psiy );
|
||||
the_output_plugin->write_dm_position(2, Psiz );
|
||||
|
||||
the_output_plugin->write_dm_velocity(0, Vx );
|
||||
the_output_plugin->write_dm_velocity(1, Vy );
|
||||
the_output_plugin->write_dm_velocity(2, Vz );
|
||||
|
||||
the_output_plugin->finalize();
|
||||
delete the_output_plugin;
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
|
60
src/output_plugin.cc
Normal file
60
src/output_plugin.cc
Normal file
|
@ -0,0 +1,60 @@
|
|||
/*
|
||||
|
||||
output.cc - This file is part of MUSIC -
|
||||
a code to generate multi-scale initial conditions
|
||||
for cosmological simulations
|
||||
|
||||
Copyright (C) 2010 Oliver Hahn
|
||||
|
||||
*/
|
||||
|
||||
#include "output_plugin.hh"
|
||||
|
||||
|
||||
std::map< std::string, output_plugin_creator *>&
|
||||
get_output_plugin_map()
|
||||
{
|
||||
static std::map< std::string, output_plugin_creator* > output_plugin_map;
|
||||
return output_plugin_map;
|
||||
}
|
||||
|
||||
void print_output_plugins()
|
||||
{
|
||||
std::map< std::string, output_plugin_creator *>& m = get_output_plugin_map();
|
||||
|
||||
std::map< std::string, output_plugin_creator *>::iterator it;
|
||||
it = m.begin();
|
||||
std::cout << " - Available output plug-ins:\n";
|
||||
while( it!=m.end() )
|
||||
{
|
||||
if( (*it).second )
|
||||
std::cout << "\t\'" << (*it).first << "\'\n";
|
||||
++it;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
output_plugin *select_output_plugin( ConfigFile& cf )
|
||||
{
|
||||
std::string formatname = cf.GetValue<std::string>( "output", "format" );
|
||||
|
||||
output_plugin_creator *the_output_plugin_creator
|
||||
= get_output_plugin_map()[ formatname ];
|
||||
|
||||
if( !the_output_plugin_creator )
|
||||
{
|
||||
std::cerr << " - Error: output plug-in \'" << formatname << "\' not found." << std::endl;
|
||||
print_output_plugins();
|
||||
throw std::runtime_error("Unknown output plug-in");
|
||||
|
||||
}else
|
||||
std::cout << " - Selecting output plug-in \'" << formatname << "\'..." << std::endl;
|
||||
|
||||
output_plugin *the_output_plugin
|
||||
= the_output_plugin_creator->create( cf );
|
||||
|
||||
return the_output_plugin;
|
||||
}
|
||||
|
||||
|
||||
|
1407
src/plugins/output_gadget2.cc
Normal file
1407
src/plugins/output_gadget2.cc
Normal file
File diff suppressed because it is too large
Load diff
Loading…
Reference in a new issue