core_app.cpp File Reference

implementation for common functions for all types of approximations More...

#include <algorithm>
#include "core_app.h"
#include "core_mesh.h"
#include "CBRNG_Random.h"

Include dependency graph for core_app.cpp:

Go to the source code of this file.

Macros
#define	ORDER   1
#define	N_MAX   1

Functions
template<typename T >
static T	mean (const std::vector< T > &data)
static double	mean (const Mesh &data)
template<typename T >
static T	std_dev (const std::vector< T > &data, T mean)
static double	std_dev (const Mesh &data, double mean)
template<typename T >
static T	min (const std::vector< T > &data)
static double	min (const Mesh &data)
template<typename T >
static T	max (const std::vector< T > &data)
static double	max (const Mesh &data)
static void	set_unpert_pos_one_par (Vec_3D< size_t > &unpert_pos, const size_t par_index, const size_t par_per_dim, const size_t Ng)
static void	set_velocity_one_par (const Vec_3D< size_t > &unpert_pos, Vec_3D< double > &displ_field, const std::vector< Mesh > &vel_field)
void	set_unpert_pos (const Sim_Param &sim, std::vector< Particle_x< double >> &particles)
void	set_unpert_pos_w_vel (const Sim_Param &sim, std::vector< Particle_v< double >> &particles, const std::vector< Mesh > &vel_field)
void	set_pert_pos (const Sim_Param &sim, const double a, std::vector< Particle_x< double >> &particles, const std::vector< Mesh > &vel_field)
void	set_pert_pos (const Sim_Param &sim, const double a, std::vector< Particle_v< double >> &particles, const std::vector< Mesh > &vel_field)
static void	gen_gauss_white_noise (const Sim_Param &sim, Mesh &rho)
static double	truncation_fce (double k, double k2_G)
static void	gen_rho_w_pow_k (const Sim_Param &sim, Mesh &rho)
void	gen_rho_dist_k (const Sim_Param &sim, Mesh &rho, const FFTW_PLAN_TYPE &p_F)
	Generate density distributions in k-space. More...
template<class T >
void	get_rho_from_par (const std::vector< T > &particles, Mesh &rho, const Sim_Param &sim)
bool	get_vel_from_par (const std::vector< Particle_v< double >> &particles, std::vector< Mesh > &vel_field, const Sim_Param &sim)
bool	get_vel_from_par (const std::vector< Particle_x< double >> &particles, std::vector< Mesh > &vel_field, const Sim_Param &sim)
void	pwr_spec_k (const Mesh &rho_k, Mesh &power_aux)
void	pwr_spec_k_init (const Mesh &rho_k, Mesh &power_aux)
void	vel_pwr_spec_k (const std::vector< Mesh > &vel_field, Mesh &power_aux)
void	gen_cqty_binned (const double x_min, const double x_max, const size_t bins_per_decade, const Mesh &qty_mesh, const size_t half_length, Data_Vec< double, 2 > &qty_binned, const double mod_q, const double mod_x)
void	gen_cqty_binned (const double x_min, const double x_max, const size_t bins_per_decade, const Mesh &qty_mesh, Data_Vec< double, 2 > &qty_binned, const double mod_q, const double mod_x)
void	gen_pow_spec_binned (const Sim_Param &sim, const Mesh &power_aux, Data_Vec< double, 2 > &pwr_spec_binned)
void	gen_pow_spec_binned_init (const Sim_Param &sim, const Mesh &power_aux, const size_t half_length, Data_Vec< double, 2 > &pwr_spec_binned)
template<class P , typename T , size_t N>
void	gen_pow_spec_binned_from_extrap (const Sim_Param &sim, const P &P_k, Data_Vec< T, N > &pwr_spec_binned)
void	gen_pot_k (const Mesh &rho_k, Mesh &pot_k)
void	gen_pot_k (Mesh &rho_k)
static double	S2_shape (const double k2, const double a)
static double	CIC_opt (Vec_3D< double > k_vec, const double a)
void	gen_displ_k_S2 (std::vector< Mesh > &vel_field, const Mesh &pot_k, const double a)
void	gen_displ_k (std::vector< Mesh > &vel_field, const Mesh &pot_k)
void	gen_displ_k_cic (std::vector< Mesh > &vel_field, const Mesh &pot_k)
void	gen_dens_binned (const Mesh &rho, std::vector< size_t > &dens_binned, const Sim_Param &sim)
template void	get_rho_from_par (const std::vector< Particle_x< double >> &, Mesh &, const Sim_Param &)
template void	get_rho_from_par (const std::vector< Particle_v< double >> &, Mesh &, const Sim_Param &)
template void	gen_pow_spec_binned_from_extrap (const Sim_Param &, const Extrap_Pk< double, 2 > &, Data_Vec< double, 2 > &)

Detailed Description

implementation for common functions for all types of approximations

Author

Michal Vrastil

Date

2018-07-11

Definition in file core_app.cpp.

Macro Definition Documentation

#define N_MAX   1

Referenced by CIC_opt().

#define ORDER   1

Definition at line 15 of file core_app.cpp.

Referenced by pwr_spec_k(), and vel_pwr_spec_k().

Function Documentation

static double CIC_opt ( Vec_3D< double >  k_vec, const double  a  )

static

Definition at line 531 of file core_app.cpp.

References N_MAX, PI, pow2(), and S2_shape().

Referenced by gen_displ_k_S2().

{
#define N_MAX 1
#ifndef N_MAX
    FTYPE_t s2 = pow2(S2_shape(k_vec.norm2(), a));
    for(unsigned int j=0; j<3; j++)
    {
        if (k_vec[j] != 0) s2 /= pow2(sin(k_vec[j] / 2) / (k_vec[j] / 2)); //W (k) for CIC (order 1)
    }
    return s2;
#else
    FTYPE_t k_n[3];
    FTYPE_t U2, U_n, G_n, k2n;
    
    G_n = 0;
    U2 = 1;
    for (unsigned int j = 0; j < 3; j++) U2 *= (1+2*pow2(cos(k_vec[j]/2)))/3; // inf sum of U_n^2 for CIC
    for (int n1 = -N_MAX; n1 < N_MAX + 1; n1++)
    {
        k_n[0] = k_vec[0] + 2 * PI*n1;
        for (int n2 = -N_MAX; n2 < N_MAX + 1; n2++)
        {
            k_n[1] = k_vec[1] + 2 * PI*n2;
            for (int n3 = -N_MAX; n3 < N_MAX + 1; n3++)
            {
                k_n[2] = k_vec[2] + 2 * PI*n3;
                U_n = 1.;
                k2n = 0;
                for(unsigned int j=0; j<3; j++)
                {
                    if (k_n[j] != 0) U_n *= sin(k_n[j] / 2) / (k_n[j] / 2);
                    k2n += pow2(k_n[j]);
                }
                U_n = pow2(U_n); // W(k) for CIC (order 1)
                if (k2n != 0)
                {
                    for(unsigned int j=0; j<3; j++)
                    {                                       
                        G_n += k_vec[j]* // i.D(k)
                        k_n[j]/k2n*pow2(S2_shape(k2n, a))* // R*(k_n) /i
                        U_n; // U_n
                    }
                }
            }
        }
    }
    if ((G_n != G_n) || (U2 != U2))
    {
        BOOST_LOG_TRIVIAL(warning) << "Gn = " << G_n << "\tU2 = " << U2 << ", k = (" << k_vec[0] << ", " << k_vec[1] << ", " << k_vec[2] << ")";
        return 1.;
    }
    return G_n/U2;
#endif
}

void gen_cqty_binned	(	const double	x_min,
		const double	x_max,
		const size_t	bins_per_decade,
		const Mesh &	qty_mesh,
		const size_t	half_length,
		Data_Vec< double, 2 > &	qty_binned,
		const double	mod_q,
		const double	mod_x
	)

Definition at line 401 of file core_app.cpp.

References Catch::Matchers::Vector::Detail::count(), Data_Vec< T, N >::erase(), Data_Vec< T, N >::fill(), Data_Vec< T, N >::resize(), Data_Vec< T, N >::size(), and x.

Referenced by gen_cqty_binned(), gen_pow_spec_binned(), and gen_pow_spec_binned_init().

{
    /* bin some complex quantity on mesh in logarithmic bins, assuming:
       Q(x) = mod_q*qty_mesh[2*i]
       x = mod_x*qty_mesh[2*i+1]
       [mesh[2*i+1]] = [x_min] = [x_max]

       return binned data in qty_binned {x, <Q(x)>}

       Note: passing length of the array for case when my mesh is bigger than data stored in there
             overloaded function exists when this is not the case
    */

    size_t req_size = (size_t)ceil(bins_per_decade*log10(x_max/x_min));
    qty_binned.resize(req_size);
    qty_binned.fill(0);
    std::vector<size_t> tmp(req_size, 0); // for counts in bins

    FTYPE_t x;
    size_t bin;
    
    /* compute sum x, Q(x), Q^2(x) in bins */
    #pragma omp parallel for private(x, bin)
    for (size_t i = 0; i < half_length; i++){
        x = qty_mesh[2*i+1];
        if ((x <x_max) && (x>=x_min)){
            bin = (size_t)((log10(x) - log10(x_min)) * bins_per_decade);
            #pragma omp atomic
            qty_binned[0][bin] += x;
            #pragma omp atomic
            qty_binned[1][bin] += qty_mesh[2*i];
            #pragma omp atomic
            tmp[bin]++;
        }
    }

    /* compute average x, Q(x) in bins */
    size_t count;
    for (size_t j = 0; j < qty_binned.size(); ){
        count = tmp[j];
        if (count){
            qty_binned[0][j] *= mod_x / count;
            qty_binned[1][j] *= mod_q / count;
            j++;
        }else{
            qty_binned.erase(j);
            tmp.erase(tmp.begin() + j);
        }
    }
}

void gen_cqty_binned	(	const double	x_min,
		const double	x_max,
		const size_t	bins_per_decade,
		const Mesh &	qty_mesh,
		Data_Vec< double, 2 > &	qty_binned,
		const double	mod_q,
		const double	mod_x
	)

Definition at line 454 of file core_app.cpp.

References gen_cqty_binned(), and Mesh_base< T >::length.

{
    gen_cqty_binned(x_min, x_max, bins_per_decade, qty_mesh, qty_mesh.length / 2, qty_binned, mod_q, mod_x);
}

void gen_dens_binned	(	const Mesh &	rho,
		std::vector< size_t > &	dens_binned,
		const Sim_Param &	sim
	)

Definition at line 627 of file core_app.cpp.

References Sim_Param::box_opt, gen_pow_spec_binned_from_extrap(), get_rho_from_par(), Mesh::N, Box_Opt::Ng_pwr, and pow().

Referenced by App_Var< T >::Impl< T >::print_density().

{
    BOOST_LOG_TRIVIAL(debug) << "Computing binned density field...";
    size_t bin;
    FTYPE_t rho_avg;
    const size_t Ng_pwr = sim.box_opt.Ng_pwr;
    const size_t N = rho.N;

    dens_binned.assign(dens_binned.size(), 0);
    
    #pragma omp parallel for private(bin, rho_avg)
    for (size_t i = 0; i < N; i+=Ng_pwr)
    {
        for (size_t j = 0; j < N; j+=Ng_pwr)
        {
            for (size_t k = 0; k < N; k+=Ng_pwr)
            {
                // Need to go through all mesh cells [i, i+Ng-1]*[j, j+Ng-1], [k, k+Ng, -1]
                rho_avg = 0;
                for (size_t ii = i; ii < i+Ng_pwr; ii++)
                {
                    for (size_t jj = j; jj  < j+Ng_pwr; jj++)
                    {
                        for (size_t kk = k; kk < k+Ng_pwr; kk++)
                        {
                            rho_avg+=rho(ii, jj, kk);
                        }
                    }
                }
                rho_avg /= pow((FTYPE_t)Ng_pwr, 3);
                bin = (size_t)((rho_avg+1)/0.1);
                if (bin >= dens_binned.size()) bin = dens_binned.size() - 1;
                // if (bin >= dens_binned.capacity()) dens_binned.resize(bin+1);
                #pragma omp atomic
                dens_binned[bin]++;
                //dens_binned[bin] += pow(sim.box_opt.Ng_pwr, 3);
            }
        }
    }
}

void gen_displ_k	(	std::vector< Mesh > &	vel_field,
		const Mesh &	pot_k
	)

Definition at line 623 of file core_app.cpp.

References gen_displ_k_S2().

Referenced by App_Var_AA::AAImpl::aa_convolution(), and App_Var< T >::Impl< T >::set_init_cond().

{gen_displ_k_S2(vel_field, pot_k, -1);}

void gen_displ_k_cic	(	std::vector< Mesh > &	vel_field,
		const Mesh &	pot_k
	)

Definition at line 625 of file core_app.cpp.

References gen_displ_k_S2().

Referenced by App_Var_Chi::ChiImpl::get_chi_force(), and App_Var< T >::pot_corr().

{gen_displ_k_S2(vel_field, pot_k, 0.);}

void gen_displ_k_S2	(	std::vector< Mesh > &	vel_field,
		const Mesh &	pot_k,
		const double	a
	)

Definition at line 586 of file core_app.cpp.

References CIC_opt(), get_k_vec(), and PI.

Referenced by gen_displ_k(), gen_displ_k_cic(), and App_Var_FP_mod::pot_corr().

{   /*
    pot_k can be Mesh of differen (bigger) size than each vel_field,
    !!!> ALL physical FACTORS ARE therefore TAKEN FROM vel_field[0] <!!!
    */
    if (a == -1) BOOST_LOG_TRIVIAL(debug) << "Computing displacement in k-space...";
    else if (a == 0) BOOST_LOG_TRIVIAL(debug) << "Computing displacement in k-space with CIC opt...";
    else BOOST_LOG_TRIVIAL(debug) << "Computing force in k-space for S2 shaped particles with CIC opt...";

    FTYPE_t opt;
    Vec_3D<int> k_vec;
    Vec_3D<FTYPE_t> k_vec_phys;
    FTYPE_t potential_tmp[2];
    
    const size_t N = vel_field[0].N; // for case when pot_k is different mesh than vel_field
    const FTYPE_t d_k = 2*PI/N;  // 2*PI/N comes from derivative WITH RESPECT to the mesh coordinates
    const size_t l_half = vel_field[0].length/2;
    
    #pragma omp parallel for private(opt, k_vec, k_vec_phys, potential_tmp)
    for(size_t i=0; i < l_half;i++)
    {
        potential_tmp[0] = pot_k[2*i]; // prevent overwriting if vel_field[0] == pot_k
        potential_tmp[1] = pot_k[2*i+1]; // prevent overwriting if vel_field[0] == pot_k
        get_k_vec(N, i, k_vec);
        k_vec_phys = d_k*k_vec; 
        // no optimalization
        if (a == -1) opt = 1.;
        // optimalization for CIC and S2 shaped particle
        else opt = CIC_opt(k_vec_phys, a);
        for(size_t j=0; j<3;j++)
        {
            vel_field[j][2*i] = k_vec_phys[j]*potential_tmp[1]*opt;
            vel_field[j][2*i+1] = -k_vec_phys[j]*potential_tmp[0]*opt;
        }
    }
}

static void gen_gauss_white_noise ( const Sim_Param &  sim, Mesh &  rho  )

static

Definition at line 165 of file core_app.cpp.

References max(), mean(), min(), Mesh::N, Mesh_base< T >::N1, Sim_Param::run_opt, Run_Opt::seed, sqrt(), and std_dev().

Referenced by gen_rho_dist_k().

{
    // Get keys for each slab in the x axis that this rank contains
    std::vector<size_t> slab_keys;
    slab_keys.resize(rho.N1);
    GetSlabKeys(slab_keys.data(), 0, rho.N1, sim.run_opt.seed);
    
    size_t ikey, index;
    FTYPE_t rn1, rn2, rn, tmp;
    const size_t N = rho.N;
        
    #pragma omp parallel for private(ikey, index, rn1, rn2, rn, tmp)
    for(size_t i=0; i < N; ++i)
    {
        ikey = slab_keys[i];
        for(size_t j=0; j < N; ++j) 
        {
            #ifndef NOISE_HALF
            for(size_t k=0; k< N; ++k)
            {
                index = j*N + k; 
                GetRandomDoublesWhiteNoise(rn1, rn2, rn, ikey, index);
                tmp = sqrt(-2*log(rn)/rn);
                rho(i, j, k) = rn2 * tmp;
            }
            #else
            for(size_t k=0; k < N/2; ++k) // go over half, use both random numbers
            {
                index = j*N + k;
                GetRandomDoublesWhiteNoise(rn1, rn2, rn, ikey, index);
                tmp = sqrt(-2*log(rn)/rn);
                rho(i, j, 2*k) = rn2 * tmp;
                rho(i, j, 2*k+1) = rn1 * tmp;
            }
            #endif
        }
    }     
    FTYPE_t t_mean;
    #ifdef CORR
    t_mean = mean(rho);
    FTYPE_t t_std_dev = std_dev(rho, t_mean);
    BOOST_LOG_TRIVIAL(debug) << "\t[mean = " << t_mean << ", stdDev = " << t_std_dev << "]\t-->";
    rho-=t_mean;
    rho/=t_std_dev;
    #endif
    
    t_mean = mean(rho);
    BOOST_LOG_TRIVIAL(debug) << "\t[mean = " << t_mean << ", stdDev = " << std_dev(rho, t_mean) << "]\t<--";
    BOOST_LOG_TRIVIAL(debug) << "\t[min = " << min(rho) << ", max = " << max(rho) << "]";
}

void gen_pot_k	(	const Mesh &	rho_k,
		Mesh &	pot_k
	)

Definition at line 496 of file core_app.cpp.

References get_k_sq(), Mesh_base< T >::length, Mesh::N, PI, and pow2().

Referenced by gen_pot_k(), anonymous_namespace{test_chameleon.cpp}::get_grav_pot(), App_Var< T >::Impl< T >::set_init_cond(), and App_Var_PM::upd_pos().

{   /*
    pot_k can be Mesh of differen (bigger) size rho_k,
    !!!> ALL physical FACTORS ARE therefore TAKEN FROM rho_k <!!!
    */
    BOOST_LOG_TRIVIAL(debug) << "Computing potential in k-space...";
    FTYPE_t k2;
    const size_t N = rho_k.N; // for case when pot_k is different mesh than vel_field
    const FTYPE_t d2_k = pow2(2*PI/N); // factor from second derivative with respect to the mesh coordinates
    const size_t l_half = rho_k.length/2;

    #pragma omp parallel for private(k2)
    for(size_t i=0; i < l_half;i++){                
        k2 = get_k_sq(N, i);
        if (k2 == 0){
            pot_k[2*i] = 0;
            pot_k[2*i+1] = 0;
        } else{
            pot_k[2*i] = -rho_k[2*i]/(k2*d2_k);
            pot_k[2*i+1] = -rho_k[2*i+1]/(k2*d2_k);
        }
    }
}

void gen_pot_k

(

Mesh &

rho_k

)

Definition at line 520 of file core_app.cpp.

References gen_pot_k().

{ gen_pot_k(rho_k, rho_k); }

void gen_pow_spec_binned	(	const Sim_Param &	sim,
		const Mesh &	power_aux,
		Data_Vec< double, 2 > &	pwr_spec_binned
	)

Definition at line 460 of file core_app.cpp.

References Out_Opt::bins_per_decade, Sim_Param::box_opt, Box_Opt::box_size, gen_cqty_binned(), Box_Opt::mesh_num_pwr, Sim_Param::out_opt, PI, and pow().

Referenced by App_Var_Chi::ChiImpl::gen_pow_spec_binned(), App_Var< T >::Impl< T >::get_binned_power_spec(), and App_Var< T >::Impl< T >::print_vel_pwr().

{
    const FTYPE_t mod_pk = pow(sim.box_opt.box_size, 3); // P(k) -> dimensionFULL!
    const FTYPE_t mod_k = 2*PI/sim.box_opt.box_size;
    BOOST_LOG_TRIVIAL(debug) << "Computing binned power spectrum...";
    gen_cqty_binned(1, sim.box_opt.mesh_num_pwr,  sim.out_opt.bins_per_decade, power_aux, pwr_spec_binned, mod_pk, mod_k);
}

template<class P , typename T , size_t N>

void gen_pow_spec_binned_from_extrap	(	const Sim_Param &	sim,
		const P &	P_k,
		Data_Vec< T, N > &	pwr_spec_binned
	)

Definition at line 478 of file core_app.cpp.

References Out_Opt::bins_per_decade, Other_par::k_print, Range::lower, Sim_Param::other_par, Sim_Param::out_opt, pow(), Data_Vec< T, N >::resize(), Data_Vec< T, N >::size(), growth_allz::T, and Range::upper.

Referenced by gen_dens_binned(), and App_Var< T >::Impl< T >::print_extrap_pwr().

{
    const T k_max = sim.other_par.k_print.upper;
    const T k_min = sim.other_par.k_print.lower;
    const T log_bin = T(1) / sim.out_opt.bins_per_decade;
    T k;
    size_t req_size = (size_t)ceil( sim.out_opt.bins_per_decade*log10(k_max/k_min));
    pwr_spec_binned.resize(req_size);

    #pragma omp parallel for private(k)
    for (size_t j = 0; j < pwr_spec_binned.size(); j++){
        k = k_min*pow(T(10), j*log_bin);
        pwr_spec_binned[0][j] = k;
        pwr_spec_binned[1][j] = P_k(k);
    }
}

template void gen_pow_spec_binned_from_extrap	(	const Sim_Param &	,
		const Extrap_Pk< double, 2 > &	,
		Data_Vec< double, 2 > &
	)

void gen_pow_spec_binned_init	(	const Sim_Param &	sim,
		const Mesh &	power_aux,
		const size_t	half_length,
		Data_Vec< double, 2 > &	pwr_spec_binned
	)

Definition at line 468 of file core_app.cpp.

References Out_Opt::bins_per_decade, Sim_Param::box_opt, Box_Opt::box_size, gen_cqty_binned(), Box_Opt::mesh_num, Sim_Param::out_opt, PI, and pow().

Referenced by App_Var< T >::Impl< T >::print_input_realisation().

{
    /* same as above but now  power_aux is storing only data [0...mesh_num], NOT mesh_num_pwr */
    const FTYPE_t mod_pk = pow(sim.box_opt.box_size, 3); // P(k) -> dimensionFULL!
    const FTYPE_t mod_k = 2*PI/sim.box_opt.box_size;
    BOOST_LOG_TRIVIAL(debug) << "Computing binned initial power spectrum...";
    gen_cqty_binned(1, sim.box_opt.mesh_num,  sim.out_opt.bins_per_decade, power_aux, half_length, pwr_spec_binned, mod_pk, mod_k);
}

void gen_rho_dist_k	(	const Sim_Param &	sim,
		Mesh &	rho,
		const FFTW_PLAN_TYPE &	p_F
	)

Generate density distributions in k-space.

Parameters

sim
rho
p_F

At first, a gaussian white noise (mean = 0, stdDev = 1) is generated, then it is convoluted with given power spectrum.

Definition at line 267 of file core_app.cpp.

References fftw_execute_dft_r2c(), gen_gauss_white_noise(), and gen_rho_w_pow_k().

Referenced by App_Var< T >::Impl< T >::set_init_cond().

{
    BOOST_LOG_TRIVIAL(debug) << "Generating gaussian white noise...";
    gen_gauss_white_noise(sim, rho);
    fftw_execute_dft_r2c(p_F, rho);

    BOOST_LOG_TRIVIAL(debug) << "Generating density distributions with given power spectrum...";
    gen_rho_w_pow_k(sim, rho);
}

static void gen_rho_w_pow_k ( const Sim_Param &  sim, Mesh &  rho  )

static

Definition at line 221 of file core_app.cpp.

References Sim_Param::box_opt, Box_Opt::box_size, Sim_Param::cosmo, get_k_sq(), Cosmo_Param::k2_G, Mesh_base< T >::length, lin_pow_spec(), Mesh::N, Run_Opt::phase, PI, pow(), Sim_Param::run_opt, sqrt(), Cosmo_Param::truncated_pk, and truncation_fce().

Referenced by gen_rho_dist_k().

{
    // extract const variables
    const FTYPE_t L = sim.box_opt.box_size;
    const FTYPE_t k0 = 2*PI/L;
    const int phase = sim.run_opt.phase ? 1 : -1;
    const size_t N = rho.N;
    const size_t len = rho.length / 2;
    const FTYPE_t mod = phase * pow(N / L, 3/2.); // pair sim, gaussian real -> fourier factor, dimension trans. Pk -> Pk*
    const bool truncation = sim.cosmo.truncated_pk;

    if(truncation)
    {
        const FTYPE_t k2_G = sim.cosmo.k2_G;
        #pragma omp parallel for
        for(size_t i=0; i < len; i++)
        {
            FTYPE_t k = k0*sqrt(get_k_sq(N, i));
            FTYPE_t tr = truncation_fce(k, k2_G);
            rho[2*i] *= mod*sqrt(lin_pow_spec(1, k, sim.cosmo))*tr;
            rho[2*i+1] *= mod*sqrt(lin_pow_spec(1, k, sim.cosmo))*tr;
        }
    }
    else
    {
        #pragma omp parallel for
        for(size_t i=0; i < len; i++)
        {
            FTYPE_t k = k0*sqrt(get_k_sq(N, i));
            rho[2*i] *= mod*sqrt(lin_pow_spec(1, k, sim.cosmo));
            rho[2*i+1] *= mod*sqrt(lin_pow_spec(1, k, sim.cosmo));
        }
    }
    
}

template<class T >

void get_rho_from_par	(	const std::vector< T > &	particles,
		Mesh &	rho,
		const Sim_Param &	sim
	)

Definition at line 278 of file core_app.cpp.

References Mesh_base< T >::assign(), assign_to(), Sim_Param::box_opt, m, Box_Opt::mesh_num, Mesh::N, Box_Opt::par_num, and pow().

Referenced by gen_dens_binned(), App_Var< T >::Impl< T >::print_output(), App_Var_Chi::ChiImpl::solve(), and App_Var_PM::upd_pos().

{
    BOOST_LOG_TRIVIAL(debug) << "Computing the density field from particle positions...";

    const size_t Np = sim.box_opt.par_num;
    if (particles.size() != Np){
        std::string msg = "Number of particles (" + std::to_string(particles.size()) + ") does not correspond with simulation parameters (" + std::to_string(Np) + ")!";
        throw std::range_error(msg);
    }
    const FTYPE_t m = pow((FTYPE_t)rho.N, 3) / Np;
    const FTYPE_t mesh_mod = (FTYPE_t)rho.N/sim.box_opt.mesh_num;
    

    rho.assign(-1.);
    
    #pragma omp parallel for
    for (size_t i = 0; i < Np; i++)
    {
        assign_to(rho, particles[i].position*mesh_mod, m);
    }
}

template void get_rho_from_par	(	const std::vector< Particle_x< double >> &	,
		Mesh &	,
		const Sim_Param &
	)

template void get_rho_from_par	(	const std::vector< Particle_v< double >> &	,
		Mesh &	,
		const Sim_Param &
	)

bool get_vel_from_par	(	const std::vector< Particle_v< double >> &	particles,
		std::vector< Mesh > &	vel_field,
		const Sim_Param &	sim
	)

Definition at line 300 of file core_app.cpp.

References assign_to(), Sim_Param::box_opt, m, Box_Opt::mesh_num, Box_Opt::mesh_num_pwr, Box_Opt::Ng_pwr, Box_Opt::par_num, and pow().

Referenced by App_Var< T >::Impl< T >::print_output().

{
    BOOST_LOG_TRIVIAL(debug) << "Computing the velocity field from particle positions...";
    const FTYPE_t mesh_mod = (FTYPE_t)sim.box_opt.mesh_num_pwr/sim.box_opt.mesh_num;
    const FTYPE_t m = pow((FTYPE_t)sim.box_opt.Ng_pwr, 3);
    const size_t Np = sim.box_opt.par_num;

    for(Mesh& field : vel_field){
        field.assign(0.);
    }
    #pragma omp parallel for
    for (size_t i = 0; i < Np; i++)
    {
        assign_to(vel_field, particles[i].position*mesh_mod, particles[i].velocity*(m*mesh_mod));
    }
    return true;
}

bool get_vel_from_par	(	const std::vector< Particle_x< double >> &	particles,
		std::vector< Mesh > &	vel_field,
		const Sim_Param &	sim
	)

Definition at line 318 of file core_app.cpp.

{
    BOOST_LOG_TRIVIAL(debug) << "WARNING! Trying to compute velocity divergence with particle positions only! Skipping...";
    return false;
}

template<typename T >

static T max ( const std::vector< T > &  data )

static

Definition at line 61 of file core_app.cpp.

Referenced by assign_to_grid(), ContinuitySolver< NDIM, T >::bfunc(), nlohmann::detail::binary_writer< BasicJsonType, CharType >::calc_bson_integer_size(), nlohmann::detail::binary_writer< BasicJsonType, CharType >::calc_bson_unsigned_size(), gen_gauss_white_noise(), cl_cmbl_bm.Power_Spectra::kappa_cl(), cl_cmblx_bm.Power_Spectra::kappa_cl(), max(), nlohmann::detail::binary_writer< BasicJsonType, CharType >::ubjson_prefix(), nlohmann::detail::binary_writer< BasicJsonType, CharType >::write_bson_integer(), nlohmann::detail::binary_writer< BasicJsonType, CharType >::write_bson_unsigned(), nlohmann::detail::binary_writer< BasicJsonType, CharType >::write_cbor(), nlohmann::detail::binary_writer< BasicJsonType, CharType >::write_msgpack(), and nlohmann::detail::binary_writer< BasicJsonType, CharType >::write_number_with_ubjson_prefix().

{
    return *std::max_element(data.begin(), data.end());
}

static double max ( const Mesh &  data )

static

Definition at line 66 of file core_app.cpp.

References Mesh_base< T >::data, and max().

                                    {
    return max(data.data);
}

template<typename T >

static T mean ( const std::vector< T > &  data )

static

Definition at line 19 of file core_app.cpp.

References growth_allz::T.

Referenced by gen_gauss_white_noise(), and mean().

{
    T tmp(0);
    
    #pragma omp parallel for reduction(+:tmp)
    for (auto it = data.begin(); it < data.end(); ++it) tmp += *it;
    
    return tmp / data.size();
}

static double mean ( const Mesh &  data )

static

Definition at line 29 of file core_app.cpp.

References Mesh_base< T >::data, and mean().

{
    return mean(data.data);
}

template<typename T >

static T min ( const std::vector< T > &  data )

static

Definition at line 51 of file core_app.cpp.

Referenced by gen_gauss_white_noise(), and min().

{
    return *std::min_element(data.begin(), data.end());
}

static double min ( const Mesh &  data )

static

Definition at line 56 of file core_app.cpp.

References Mesh_base< T >::data, and min().

                                    {
    return min(data.data);
}

void pwr_spec_k	(	const Mesh &	rho_k,
		Mesh &	power_aux
	)

Definition at line 324 of file core_app.cpp.

References get_k_vec(), Mesh_base< T >::length, Mesh::N, halomod_bm::NM, ORDER, PI, and pow().

Referenced by App_Var< T >::Impl< T >::get_binned_power_spec().

{
    /* Computing the power spectrum P(k)/L^3 -- dimensionLESS!

    > in real part [even] of power_aux is stored pk, in imaginary [odd] dimensionLESS k
    > preserve values in rho_k
    > as power_aux can be Mesh of different (bigger) size than rho_k, all sizes / lengths are taken from rho_k
    */
    
    FTYPE_t w_k;
    Vec_3D<int> k_vec;
    const size_t NM = rho_k.N;
    const size_t half_length = rho_k.length / 2;

    #pragma omp parallel for private(w_k, k_vec)
    for(size_t i=0; i < half_length;i++)
    {
        w_k = 1.;
        get_k_vec(NM, i, k_vec);
        for (unsigned int j = 0; j < 3; j++) if (k_vec[j] != 0) w_k *= pow(sin(k_vec[j]*PI/NM)/(k_vec[j]*PI/NM), ORDER + 1);
        power_aux[2*i] = (rho_k[2*i]*rho_k[2*i] + rho_k[2*i+1]*rho_k[2*i+1])/(w_k*w_k);
        power_aux[2*i+1] = k_vec.norm();
    }
}

void pwr_spec_k_init	(	const Mesh &	rho_k,
		Mesh &	power_aux
	)

Definition at line 349 of file core_app.cpp.

References get_k_vec(), Mesh_base< T >::length, Mesh::N, halomod_bm::NM, and pow2().

Referenced by App_Var_Chi::ChiImpl::gen_pow_spec_binned(), and App_Var< T >::Impl< T >::print_input_realisation().

{
    /* same as above but now there is NO w_k correction */

    Vec_3D<int> k_vec;
    const size_t NM = rho_k.N;
    const size_t half_length = rho_k.length / 2;

    #pragma omp parallel for private(k_vec)
    for(size_t i=0; i < half_length;i++)
    {
        get_k_vec(NM, i, k_vec);
        power_aux[2*i] = pow2(rho_k[2*i]) + pow2(rho_k[2*i+1]);
        power_aux[2*i+1] = k_vec.norm();
    }
}

static double S2_shape ( const double  k2, const double  a  )

static

Definition at line 522 of file core_app.cpp.

References pow(), and sqrt().

Referenced by CIC_opt().

{
    if (a == 0) return 1.;
    
    FTYPE_t t = sqrt(k2)*a / 2;
    if (t == 0) return 1.;
    return 12 / pow(t, 4)*(2 - 2 * cos(t) - t*sin(t));
}

void set_pert_pos	(	const Sim_Param &	sim,
		const double	a,
		std::vector< Particle_x< double >> &	particles,
		const std::vector< Mesh > &	vel_field
	)

Definition at line 114 of file core_app.cpp.

References Sim_Param::box_opt, Sim_Param::cosmo, get_per(), growth_allz::growth_factor, Box_Opt::mesh_num, Box_Opt::Ng, Box_Opt::par_num, Box_Opt::par_num_1d, set_unpert_pos_one_par(), and set_velocity_one_par().

Referenced by App_Var< T >::Impl< T >::set_init_pos(), and App_Var_ZA::upd_pos().

{
    BOOST_LOG_TRIVIAL(debug) << "Setting initial positions of particles...";
    Vec_3D<size_t> unpert_pos;
    Vec_3D<FTYPE_t> velocity;
    Vec_3D<FTYPE_t> pert_pos;
    
    const size_t par_per_dim = sim.box_opt.par_num_1d;
    const size_t Ng = sim.box_opt.Ng;
    const size_t Nm = sim.box_opt.mesh_num;
    const size_t Np = sim.box_opt.par_num;

    const FTYPE_t D = growth_factor(a, sim.cosmo); // growth factor
    
    #pragma omp parallel for private(unpert_pos, velocity, pert_pos)
    for(size_t i=0; i< Np; i++)
    {
        set_unpert_pos_one_par(unpert_pos, i, par_per_dim, Ng);
        set_velocity_one_par(unpert_pos, velocity, vel_field);
        pert_pos = velocity*D + unpert_pos;
        get_per(pert_pos, Nm);
        particles[i] = Particle_x<FTYPE_t>(pert_pos);       
    }
}

void set_pert_pos	(	const Sim_Param &	sim,
		const double	a,
		std::vector< Particle_v< double >> &	particles,
		const std::vector< Mesh > &	vel_field
	)

Definition at line 139 of file core_app.cpp.

References Sim_Param::box_opt, Sim_Param::cosmo, get_per(), growth_change(), growth_allz::growth_factor, Box_Opt::mesh_num, Box_Opt::Ng, Box_Opt::par_num, Box_Opt::par_num_1d, set_unpert_pos_one_par(), and set_velocity_one_par().

{
    BOOST_LOG_TRIVIAL(debug) << "Setting initial positions and velocitis of particles...";
    Vec_3D<size_t> unpert_pos;
    Vec_3D<FTYPE_t> velocity;
    Vec_3D<FTYPE_t> pert_pos;
    
    const size_t par_per_dim = sim.box_opt.par_num_1d;
    const size_t Ng = sim.box_opt.Ng;
    const size_t Nm = sim.box_opt.mesh_num;
    const size_t Np = sim.box_opt.par_num;

    const FTYPE_t D = growth_factor(a, sim.cosmo); // growth factor
    const FTYPE_t dDda = growth_change(a, sim.cosmo); // dD / da

    #pragma omp parallel for private(unpert_pos, velocity, pert_pos)
    for(size_t i=0; i< Np; i++)
    {
        set_unpert_pos_one_par(unpert_pos, i, par_per_dim, Ng);
        set_velocity_one_par(unpert_pos, velocity, vel_field);
        pert_pos = velocity*D + unpert_pos;
        get_per(pert_pos, Nm);
        particles[i] = Particle_v<FTYPE_t>(pert_pos, velocity*dDda);        
    }
}

void set_unpert_pos	(	const Sim_Param &	sim,
		std::vector< Particle_x< double >> &	particles
	)

Definition at line 82 of file core_app.cpp.

References Sim_Param::box_opt, Box_Opt::Ng, Box_Opt::par_num, Box_Opt::par_num_1d, and set_unpert_pos_one_par().

{
    Vec_3D<size_t> unpert_pos;
    const size_t par_per_dim = sim.box_opt.par_num_1d;
    const size_t Ng = sim.box_opt.Ng;
    const size_t Np = sim.box_opt.par_num;
    
    #pragma omp parallel for private(unpert_pos)
    for(size_t i=0; i< Np; i++)
    {
        set_unpert_pos_one_par(unpert_pos, i, par_per_dim, Ng);     
        particles[i] = Particle_x<FTYPE_t>(unpert_pos);
    }
}

static void set_unpert_pos_one_par ( Vec_3D< size_t > &  unpert_pos, const size_t  par_index, const size_t  par_per_dim, const size_t  Ng  )

static

Definition at line 70 of file core_app.cpp.

Referenced by set_pert_pos(), set_unpert_pos(), and set_unpert_pos_w_vel().

{
    unpert_pos[0] = (par_index / (par_per_dim * par_per_dim)) * Ng;
    unpert_pos[1] = ((par_index / par_per_dim) % par_per_dim) * Ng;
    unpert_pos[2] = (par_index % par_per_dim) * Ng;
}

void set_unpert_pos_w_vel	(	const Sim_Param &	sim,
		std::vector< Particle_v< double >> &	particles,
		const std::vector< Mesh > &	vel_field
	)

Definition at line 97 of file core_app.cpp.

References Sim_Param::box_opt, Box_Opt::Ng, Box_Opt::par_num, Box_Opt::par_num_1d, set_unpert_pos_one_par(), and set_velocity_one_par().

{
    Vec_3D<size_t> unpert_pos;
    Vec_3D<FTYPE_t> velocity;
    const size_t par_per_dim = sim.box_opt.par_num_1d;
    const size_t Ng = sim.box_opt.Ng;
    
    const size_t Np = sim.box_opt.par_num;
    #pragma omp parallel for private(unpert_pos, velocity)
    for(size_t i=0; i< Np; i++)
    {
        set_unpert_pos_one_par(unpert_pos, i, par_per_dim, Ng);
        set_velocity_one_par(unpert_pos, velocity, vel_field);
        particles[i] = Particle_v<FTYPE_t>(unpert_pos, velocity);
    }
}

static void set_velocity_one_par ( const Vec_3D< size_t > &  unpert_pos, Vec_3D< double > &  displ_field, const std::vector< Mesh > &  vel_field  )

static

Definition at line 77 of file core_app.cpp.

Referenced by set_pert_pos(), and set_unpert_pos_w_vel().

{
    for (size_t i = 0; i < 3; i++) displ_field[i] = vel_field[i](unpert_pos);
}

template<typename T >

static T std_dev ( const std::vector< T > &  data, T  mean  )

static

Definition at line 35 of file core_app.cpp.

References pow2(), sqrt(), and growth_allz::T.

Referenced by gen_gauss_white_noise(), and std_dev().

{
    T tmp(0);
    
    #pragma omp parallel for reduction(+:tmp)
    for (auto it = data.begin(); it < data.end(); ++it) tmp += pow2(*it-mean);
    
    return sqrt(tmp / data.size());
}

static double std_dev ( const Mesh &  data, double  mean  )

static

Definition at line 45 of file core_app.cpp.

References Mesh_base< T >::data, and std_dev().

{
    return std_dev(data.data, mean);
}

static double truncation_fce ( double  k, double  k2_G  )

static

Definition at line 216 of file core_app.cpp.

Referenced by gen_rho_w_pow_k().

{
    return exp(-k*k/k2_G);
}

void vel_pwr_spec_k	(	const std::vector< Mesh > &	vel_field,
		Mesh &	power_aux
	)

Definition at line 366 of file core_app.cpp.

References get_k_vec(), halomod_bm::NM, ORDER, PI, and pow().

Referenced by App_Var< T >::Impl< T >::print_vel_pwr().

{
    /* Computing the velocity power spectrum divergence P(k)/L^3 -- dimensionLESS!

    > in real part [even] of power_aux is stored pk, in imaginary [odd] dimensionLESS k
    > preserve values in rho_k
    > as power_aux can be Mesh of different (bigger) size than rho_k, all sizes / lengths are taken from rho_k
    */
    
    FTYPE_t w_k;
    Vec_3D<int> k_vec;

    const size_t NM = vel_field[0].N;
    const size_t half_length = vel_field[0].length / 2;

    FTYPE_t vel_div_re, vel_div_im, k; // temporary store of Pk in case vel_field[0] = power_aux

    #pragma omp parallel for private(w_k, k_vec, k, vel_div_re, vel_div_im)
    for(size_t i=0; i < half_length; i++)
    {
        w_k = 1.;
        vel_div_re = vel_div_im = 0;
        get_k_vec(NM, i, k_vec);
        for (unsigned int j = 0; j < 3; j++){
            k = k_vec[j]*2*PI / NM;
            if (k != 0) w_k *= pow(sin(k/2)/(k/2), ORDER + 1);
            vel_div_re += vel_field[j][2*i]*k; // do not care about Re <-> Im in 2*PI*i/N, norm only
            vel_div_im += vel_field[j][2*i+1]*k;
        } 
        
        power_aux[2*i] = (vel_div_re*vel_div_re + vel_div_im*vel_div_im)/(w_k*w_k);
        power_aux[2*i+1] = k_vec.norm();
    }
}