core_app.h File Reference

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

#include "stdafx.h"
#include "app_var.hpp"
#include "core_power.h"
#include "params.hpp"
#include "precision.hpp"
#include "class_particles.hpp"

Include dependency graph for core_app.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
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)
void	gen_rho_dist_k (const Sim_Param &sim, Mesh &rho, const FFTW_PLAN_TYPE &p_F)
	Generate density distributions in k-space. More...
void	gen_pot_k (const Mesh &rho_k, Mesh &pot_k)
void	gen_pot_k (Mesh &rho_k)
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_displ_k_S2 (std::vector< Mesh > &vel_field, const Mesh &pot_k, const double a)
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_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_dens_binned (const Mesh &rho, std::vector< size_t > &dens_binned, const Sim_Param &sim)

Detailed Description

interface for common functions for all types of approximations

Author

Michal Vrastil

Date

2018-07-11

Definition in file core_app.h.

Function Documentation

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;
        }
    }
}

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);
    }
}

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);
}

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);
    }
}

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;
}

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();
    }
}

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);
    }
}

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);
    }
}

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();
    }
}