2019-05-09 21:41:54 +02:00
|
|
|
#pragma once
|
2019-05-07 01:05:16 +02:00
|
|
|
|
|
|
|
#include <cmath>
|
|
|
|
#include <array>
|
|
|
|
#include <vector>
|
|
|
|
|
|
|
|
#include <vec3.hh>
|
|
|
|
#include <general.hh>
|
|
|
|
#include <bounding_box.hh>
|
2019-09-25 21:46:17 +02:00
|
|
|
#include <typeinfo>
|
2019-05-07 01:05:16 +02:00
|
|
|
|
|
|
|
enum space_t
|
|
|
|
{
|
|
|
|
kspace_id,
|
|
|
|
rspace_id
|
|
|
|
};
|
|
|
|
|
2019-05-10 04:48:35 +02:00
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
#ifdef USE_MPI
|
|
|
|
template <typename data_t, bool bdistributed=true>
|
|
|
|
#else
|
|
|
|
template <typename data_t, bool bdistributed=false>
|
|
|
|
#endif
|
2019-05-07 01:05:16 +02:00
|
|
|
class Grid_FFT
|
|
|
|
{
|
2019-05-09 21:41:54 +02:00
|
|
|
protected:
|
2019-05-07 01:05:16 +02:00
|
|
|
#if defined(USE_MPI)
|
2019-11-14 15:36:39 +01:00
|
|
|
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;
|
2019-05-07 01:05:16 +02:00
|
|
|
#endif
|
2019-11-14 15:36:39 +01:00
|
|
|
using grid_fft_t = Grid_FFT<data_t,bdistributed>;
|
2019-05-07 01:05:16 +02:00
|
|
|
public:
|
2019-05-09 21:41:54 +02:00
|
|
|
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_;
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
space_t space_;
|
|
|
|
data_t *data_;
|
|
|
|
ccomplex_t *cdata_;
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
bounding_box<size_t> global_range_;
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
fftw_plan_t plan_, iplan_;
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
real_t fft_norm_fac_;
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
ptrdiff_t local_0_start_, local_1_start_;
|
|
|
|
ptrdiff_t local_0_size_, local_1_size_;
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
Grid_FFT(const std::array<size_t, 3> &N, const std::array<real_t, 3> &L, space_t initialspace = rspace_id)
|
2019-08-05 15:04:50 +02:00
|
|
|
: n_(N), length_(L), space_(initialspace), data_(nullptr), cdata_(nullptr)
|
2019-05-09 21:41:54 +02:00
|
|
|
{
|
|
|
|
//invalidated = true;
|
|
|
|
this->Setup();
|
|
|
|
}
|
|
|
|
|
2019-05-21 00:24:09 +02:00
|
|
|
// avoid implicit copying of data
|
2019-11-14 15:36:39 +01:00
|
|
|
Grid_FFT(const grid_fft_t &g) = delete;
|
2019-05-09 21:41:54 +02:00
|
|
|
|
|
|
|
~Grid_FFT()
|
|
|
|
{
|
|
|
|
if (data_ != nullptr)
|
|
|
|
{
|
|
|
|
fftw_free(data_);
|
|
|
|
}
|
|
|
|
}
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
const grid_fft_t *get_grid(size_t ilevel) const { return this; }
|
2019-08-05 15:04:50 +02:00
|
|
|
|
2019-11-15 22:32:09 +01:00
|
|
|
bool is_distributed( void ) const { return bdistributed; }
|
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
void Setup();
|
|
|
|
|
2019-12-01 11:28:17 +01:00
|
|
|
//! return the number of data_t elements that we store in the container
|
|
|
|
size_t memsize( void ) const { return ntot_; }
|
|
|
|
|
2019-08-12 00:14:30 +02:00
|
|
|
//! return the (local) size of dimension i
|
2019-05-09 21:41:54 +02:00
|
|
|
size_t size(size_t i) const { return sizes_[i]; }
|
|
|
|
|
2019-08-12 15:35:52 +02:00
|
|
|
//! return the (global) size of dimension i
|
|
|
|
size_t global_size(size_t i) const { return n_[i]; }
|
|
|
|
|
2019-08-12 00:14:30 +02:00
|
|
|
//! return locally stored number of elements of field
|
2019-10-09 17:25:15 +02:00
|
|
|
size_t local_size(void) const { return local_0_size_ * n_[1] * n_[2]; }
|
2019-07-31 11:57:40 +02:00
|
|
|
|
2019-08-12 00:14:30 +02:00
|
|
|
//! return a bounding box of the global extent of the field
|
2019-10-09 17:25:15 +02:00
|
|
|
const bounding_box<size_t> &get_global_range(void) const
|
2019-05-20 17:23:52 +02:00
|
|
|
{
|
|
|
|
return global_range_;
|
|
|
|
}
|
|
|
|
|
2019-08-12 00:14:30 +02:00
|
|
|
//! set all field elements to zero
|
2019-05-09 21:41:54 +02:00
|
|
|
void zero()
|
|
|
|
{
|
2019-10-09 17:25:15 +02:00
|
|
|
#pragma omp parallel for
|
2019-05-09 21:41:54 +02:00
|
|
|
for (size_t i = 0; i < ntot_; ++i)
|
|
|
|
data_[i] = 0.0;
|
|
|
|
}
|
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
void copy_from(const grid_fft_t &g)
|
2019-10-09 17:25:15 +02:00
|
|
|
{
|
2019-08-05 15:04:50 +02:00
|
|
|
// make sure the two fields are in the same space
|
2019-10-09 17:25:15 +02:00
|
|
|
if (g.space_ != this->space_)
|
|
|
|
{
|
|
|
|
if (this->space_ == kspace_id)
|
|
|
|
this->FourierTransformBackward(false);
|
|
|
|
else
|
|
|
|
this->FourierTransformForward(false);
|
2019-08-05 15:04:50 +02:00
|
|
|
}
|
|
|
|
|
|
|
|
// make sure the two fields have the same dimensions
|
2019-10-09 17:25:15 +02:00
|
|
|
assert(this->n_[0] == g.n_[0]);
|
|
|
|
assert(this->n_[1] == g.n_[1]);
|
|
|
|
assert(this->n_[2] == g.n_[2]);
|
2019-08-05 15:04:50 +02:00
|
|
|
|
2019-10-09 17:25:15 +02:00
|
|
|
// now we can copy all the data over
|
|
|
|
#pragma omp parallel for
|
2019-08-05 15:04:50 +02:00
|
|
|
for (size_t i = 0; i < ntot_; ++i)
|
|
|
|
data_[i] = g.data_[i];
|
|
|
|
}
|
|
|
|
|
2019-10-09 17:25:15 +02:00
|
|
|
data_t &operator[](size_t i)
|
|
|
|
{
|
2019-05-23 14:53:11 +02:00
|
|
|
return data_[i];
|
|
|
|
}
|
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
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];
|
2019-05-07 01:05:16 +02:00
|
|
|
}
|
2019-05-09 21:41:54 +02:00
|
|
|
|
|
|
|
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>
|
2019-05-15 05:30:47 +02:00
|
|
|
vec3<ft> get_r(const size_t i, const size_t j, const size_t k) const
|
2019-05-09 21:41:54 +02:00
|
|
|
{
|
|
|
|
vec3<ft> rr;
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
rr[0] = real_t(i + local_0_start_) * dx_[0];
|
|
|
|
rr[1] = real_t(j) * dx_[1];
|
|
|
|
rr[2] = real_t(k) * dx_[2];
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
return rr;
|
|
|
|
}
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2019-07-31 11:57:40 +02:00
|
|
|
template <typename ft>
|
|
|
|
vec3<ft> get_unit_r(const size_t i, const size_t j, const size_t k) const
|
|
|
|
{
|
|
|
|
vec3<ft> rr;
|
|
|
|
|
|
|
|
rr[0] = real_t(i + local_0_start_) / real_t(n_[0]);
|
|
|
|
rr[1] = real_t(j) / real_t(n_[1]);
|
|
|
|
rr[2] = real_t(k) / real_t(n_[2]);
|
|
|
|
|
|
|
|
return rr;
|
|
|
|
}
|
2019-10-09 17:25:15 +02:00
|
|
|
|
2019-07-31 19:47:47 +02:00
|
|
|
template <typename ft>
|
2019-11-01 04:47:02 +01:00
|
|
|
vec3<ft> get_unit_r_shifted(const size_t i, const size_t j, const size_t k, const vec3<real_t> s) const
|
2019-07-31 19:47:47 +02:00
|
|
|
{
|
|
|
|
vec3<ft> rr;
|
|
|
|
|
2019-11-01 04:47:02 +01:00
|
|
|
rr[0] = (real_t(i + local_0_start_) + s.x) / real_t(n_[0]);
|
|
|
|
rr[1] = (real_t(j) + s.y) / real_t(n_[1]);
|
|
|
|
rr[2] = (real_t(k) + s.z) / real_t(n_[2]);
|
2019-07-31 19:47:47 +02:00
|
|
|
|
|
|
|
return rr;
|
|
|
|
}
|
2019-07-31 11:57:40 +02:00
|
|
|
|
2019-10-09 17:25:15 +02:00
|
|
|
vec3<size_t> get_cell_idx_3d(const size_t i, const size_t j, const size_t k) const
|
|
|
|
{
|
|
|
|
return vec3<size_t>({i + local_0_start_, j, k});
|
2019-05-09 21:41:54 +02:00
|
|
|
}
|
|
|
|
|
2019-10-09 17:25:15 +02:00
|
|
|
size_t get_cell_idx_1d(const size_t i, const size_t j, const size_t k) const
|
|
|
|
{
|
|
|
|
return ((i + local_0_start_) * size(1) + j) * size(2) + k;
|
|
|
|
}
|
|
|
|
|
|
|
|
size_t count_leaf_cells(int, int) const
|
|
|
|
{
|
|
|
|
return n_[0] * n_[1] * n_[2];
|
|
|
|
}
|
|
|
|
|
|
|
|
real_t get_dx(int idim) const
|
|
|
|
{
|
2019-07-31 19:47:47 +02:00
|
|
|
return dx_[idim];
|
|
|
|
}
|
|
|
|
|
2019-10-09 17:25:15 +02:00
|
|
|
const std::array<real_t, 3> &get_dx(void) const
|
|
|
|
{
|
2019-07-31 19:47:47 +02:00
|
|
|
return dx_;
|
|
|
|
}
|
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
template <typename ft>
|
2019-05-12 18:28:37 +02:00
|
|
|
vec3<ft> get_k(const size_t i, const size_t j, const size_t k) const
|
2019-05-09 21:41:54 +02:00
|
|
|
{
|
|
|
|
vec3<ft> kk;
|
2019-11-14 15:36:39 +01:00
|
|
|
if( bdistributed ){
|
|
|
|
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];
|
|
|
|
}
|
2019-05-09 21:41:54 +02:00
|
|
|
kk[2] = (real_t(k) - real_t(k > nhalf_[2]) * n_[2]) * kfac_[2];
|
|
|
|
|
|
|
|
return kk;
|
|
|
|
}
|
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
std::array<size_t,3> get_k3(const size_t i, const size_t j, const size_t k) const
|
|
|
|
{
|
|
|
|
return bdistributed? std::array<size_t,3>({j,i+local_1_start_,k}) : std::array<size_t,3>({i,j,k});
|
|
|
|
}
|
|
|
|
|
2019-11-05 19:14:14 +01:00
|
|
|
data_t get_cic( const vec3<real_t>& v ) const{
|
|
|
|
// warning! this doesn't work with MPI
|
|
|
|
vec3<real_t> x({std::fmod(v.x/length_[0]+1.0,1.0)*n_[0],
|
|
|
|
std::fmod(v.y/length_[1]+1.0,1.0)*n_[1],
|
|
|
|
std::fmod(v.z/length_[2]+1.0,1.0)*n_[2] });
|
|
|
|
size_t ix = static_cast<size_t>(x.x);
|
|
|
|
size_t iy = static_cast<size_t>(x.y);
|
|
|
|
size_t iz = static_cast<size_t>(x.z);
|
|
|
|
real_t dx = x.x-real_t(ix), tx = 1.0-dx;
|
|
|
|
real_t dy = x.y-real_t(iy), ty = 1.0-dy;
|
|
|
|
real_t dz = x.z-real_t(iz), tz = 1.0-dz;
|
|
|
|
size_t ix1 = (ix+1)%n_[0];
|
|
|
|
size_t iy1 = (iy+1)%n_[1];
|
|
|
|
size_t iz1 = (iz+1)%n_[2];
|
|
|
|
data_t val = 0.0;
|
|
|
|
val += this->relem(ix ,iy ,iz ) * tx * ty * tz;
|
|
|
|
val += this->relem(ix ,iy ,iz1) * tx * ty * dz;
|
|
|
|
val += this->relem(ix ,iy1,iz ) * tx * dy * tz;
|
|
|
|
val += this->relem(ix ,iy1,iz1) * tx * dy * dz;
|
|
|
|
val += this->relem(ix1,iy ,iz ) * dx * ty * tz;
|
|
|
|
val += this->relem(ix1,iy ,iz1) * dx * ty * dz;
|
|
|
|
val += this->relem(ix1,iy1,iz ) * dx * dy * tz;
|
|
|
|
val += this->relem(ix1,iy1,iz1) * dx * dy * dz;
|
|
|
|
return val;
|
|
|
|
}
|
|
|
|
|
|
|
|
ccomplex_t get_cic_kspace( const vec3<real_t>& x ) const{
|
|
|
|
// warning! this doesn't work with MPI
|
|
|
|
size_t ix = static_cast<size_t>(x.x);
|
|
|
|
size_t iy = static_cast<size_t>(x.y);
|
|
|
|
size_t iz = std::min(static_cast<size_t>(x.z),size(2)-1); //static_cast<size_t>(x.z);
|
|
|
|
real_t dx = x.x-real_t(ix), tx = 1.0-dx;
|
|
|
|
real_t dy = x.y-real_t(iy), ty = 1.0-dy;
|
|
|
|
real_t dz = x.z-real_t(iz), tz = 1.0-dz;
|
|
|
|
size_t ix1 = (ix+1)%size(0);
|
|
|
|
size_t iy1 = (iy+1)%size(1);
|
|
|
|
size_t iz1 = std::min((iz+1),size(2)-1);
|
|
|
|
ccomplex_t val = 0.0;
|
|
|
|
val += this->kelem(ix ,iy ,iz ) * tx * ty * tz;
|
|
|
|
val += this->kelem(ix ,iy ,iz1) * tx * ty * dz;
|
|
|
|
val += this->kelem(ix ,iy1,iz ) * tx * dy * tz;
|
|
|
|
val += this->kelem(ix ,iy1,iz1) * tx * dy * dz;
|
|
|
|
val += this->kelem(ix1,iy ,iz ) * dx * ty * tz;
|
|
|
|
val += this->kelem(ix1,iy ,iz1) * dx * ty * dz;
|
|
|
|
val += this->kelem(ix1,iy1,iz ) * dx * dy * tz;
|
|
|
|
val += this->kelem(ix1,iy1,iz1) * dx * dy * dz;
|
|
|
|
// if( val != val ){
|
|
|
|
//auto k = this->get_k<real_t>(ix,iy,iz);
|
|
|
|
//std::cerr << ix << " " << iy << " " << iz << " " << val << " " << this->gradient(0,{ix,iy,iz}) << " " << this->gradient(1,{ix,iy,iz}) << " " << this->gradient(2,{ix,iy,iz}) << std::endl;
|
|
|
|
// }
|
|
|
|
return val;
|
|
|
|
}
|
|
|
|
|
2019-10-15 19:48:38 +02:00
|
|
|
inline ccomplex_t gradient( const int idim, std::array<size_t,3> ijk ) const
|
|
|
|
{
|
2019-11-14 15:36:39 +01:00
|
|
|
if( bdistributed ){
|
|
|
|
ijk[0] += local_1_start_;
|
|
|
|
std::swap(ijk[0],ijk[1]);
|
|
|
|
}
|
2019-10-15 19:48:38 +02:00
|
|
|
real_t rgrad =
|
|
|
|
(ijk[idim]!=nhalf_[idim])? (real_t(ijk[idim]) - real_t(ijk[idim] > nhalf_[idim]) * n_[idim]) * kfac_[idim] : 0.0;
|
|
|
|
return ccomplex_t(0.0,rgrad);
|
|
|
|
}
|
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
grid_fft_t &operator*=(data_t x)
|
2019-10-09 17:25:15 +02:00
|
|
|
{
|
|
|
|
if (space_ == kspace_id)
|
|
|
|
{
|
|
|
|
this->apply_function_k([&](ccomplex_t &f) { return f * x; });
|
|
|
|
}
|
|
|
|
else
|
|
|
|
{
|
|
|
|
this->apply_function_r([&](data_t &f) { return f * x; });
|
2019-05-19 10:02:32 +02:00
|
|
|
}
|
|
|
|
return *this;
|
|
|
|
}
|
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
grid_fft_t &operator/=(data_t x)
|
2019-10-09 17:25:15 +02:00
|
|
|
{
|
|
|
|
if (space_ == kspace_id)
|
|
|
|
{
|
|
|
|
this->apply_function_k([&](ccomplex_t &f) { return f / x; });
|
|
|
|
}
|
|
|
|
else
|
|
|
|
{
|
|
|
|
this->apply_function_r([&](data_t &f) { return f / x; });
|
2019-05-19 10:02:32 +02:00
|
|
|
}
|
|
|
|
return *this;
|
|
|
|
}
|
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
grid_fft_t &apply_Laplacian(void)
|
2019-10-09 17:25:15 +02:00
|
|
|
{
|
2019-05-19 12:04:46 +02:00
|
|
|
this->FourierTransformForward();
|
|
|
|
this->apply_function_k_dep([&](auto x, auto k) {
|
|
|
|
real_t kmod2 = k.norm_squared();
|
2019-10-09 17:25:15 +02:00
|
|
|
return -x * kmod2;
|
2019-05-19 12:04:46 +02:00
|
|
|
});
|
|
|
|
this->zero_DC_mode();
|
|
|
|
return *this;
|
|
|
|
}
|
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
grid_fft_t &apply_negative_Laplacian(void)
|
2019-10-09 17:25:15 +02:00
|
|
|
{
|
2019-08-08 16:03:30 +02:00
|
|
|
this->FourierTransformForward();
|
|
|
|
this->apply_function_k_dep([&](auto x, auto k) {
|
|
|
|
real_t kmod2 = k.norm_squared();
|
2019-10-09 17:25:15 +02:00
|
|
|
return x * kmod2;
|
2019-08-08 16:03:30 +02:00
|
|
|
});
|
|
|
|
this->zero_DC_mode();
|
|
|
|
return *this;
|
|
|
|
}
|
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
grid_fft_t &apply_InverseLaplacian(void)
|
2019-10-09 17:25:15 +02:00
|
|
|
{
|
2019-05-19 12:04:46 +02:00
|
|
|
this->FourierTransformForward();
|
|
|
|
this->apply_function_k_dep([&](auto x, auto k) {
|
|
|
|
real_t kmod2 = k.norm_squared();
|
2019-10-09 17:25:15 +02:00
|
|
|
return -x / kmod2;
|
2019-05-19 12:04:46 +02:00
|
|
|
});
|
|
|
|
this->zero_DC_mode();
|
|
|
|
return *this;
|
|
|
|
}
|
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
template <typename functional>
|
|
|
|
void apply_function_k(const functional &f)
|
2019-05-07 01:05:16 +02:00
|
|
|
{
|
2019-05-09 21:41:54 +02:00
|
|
|
#pragma omp parallel for
|
|
|
|
for (size_t i = 0; i < sizes_[0]; ++i)
|
2019-05-07 01:05:16 +02:00
|
|
|
{
|
2019-05-09 21:41:54 +02:00
|
|
|
for (size_t j = 0; j < sizes_[1]; ++j)
|
2019-05-07 01:05:16 +02:00
|
|
|
{
|
2019-05-09 21:41:54 +02:00
|
|
|
for (size_t k = 0; k < sizes_[2]; ++k)
|
|
|
|
{
|
|
|
|
auto &elem = this->kelem(i, j, k);
|
|
|
|
elem = f(elem);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
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);
|
|
|
|
}
|
2019-05-07 01:05:16 +02:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
double compute_2norm(void)
|
|
|
|
{
|
|
|
|
real_t sum1{0.0};
|
2019-11-14 15:36:39 +01:00
|
|
|
#pragma omp parallel for reduction(+ : sum1)
|
2019-05-09 21:41:54 +02:00
|
|
|
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];
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
return sum1;
|
|
|
|
}
|
|
|
|
|
|
|
|
double std(void)
|
|
|
|
{
|
2019-08-07 20:16:50 +02:00
|
|
|
double sum1{0.0}, sum2{0.0};
|
|
|
|
size_t count{0};
|
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
#pragma omp parallel for reduction(+ : sum1, sum2)
|
2019-05-09 21:41:54 +02:00
|
|
|
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;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2019-08-07 20:16:50 +02:00
|
|
|
count = sizes_[0] * sizes_[1] * sizes_[2];
|
2019-05-09 21:41:54 +02:00
|
|
|
|
2019-08-07 20:16:50 +02:00
|
|
|
#ifdef USE_MPI
|
2019-11-14 15:36:39 +01:00
|
|
|
if( bdistributed ){
|
|
|
|
double globsum1{0.0}, globsum2{0.0};
|
|
|
|
size_t globcount{0};
|
2019-08-07 20:16:50 +02:00
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
MPI_Allreduce(reinterpret_cast<const void *>(&sum1),
|
|
|
|
reinterpret_cast<void *>(&globsum1),
|
|
|
|
1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
|
2019-08-07 20:16:50 +02:00
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
MPI_Allreduce(reinterpret_cast<const void *>(&sum2),
|
|
|
|
reinterpret_cast<void *>(&globsum2),
|
|
|
|
1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
|
2019-08-07 20:16:50 +02:00
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
MPI_Allreduce(reinterpret_cast<const void *>(&count),
|
|
|
|
reinterpret_cast<void *>(&globcount),
|
|
|
|
1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, MPI_COMM_WORLD);
|
2019-08-07 20:16:50 +02:00
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
sum1 = globsum1;
|
|
|
|
sum2 = globsum2;
|
|
|
|
count = globcount;
|
|
|
|
}
|
2019-10-09 17:25:15 +02:00
|
|
|
#endif
|
2019-08-07 20:16:50 +02:00
|
|
|
sum1 /= count;
|
|
|
|
sum2 /= count;
|
2019-05-09 21:41:54 +02:00
|
|
|
|
|
|
|
return std::sqrt(sum2 - sum1 * sum1);
|
|
|
|
}
|
2019-08-07 20:16:50 +02:00
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
double mean(void)
|
|
|
|
{
|
2019-08-07 20:16:50 +02:00
|
|
|
double sum1{0.0};
|
|
|
|
size_t count{0};
|
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
#pragma omp parallel for reduction(+ : sum1)
|
2019-05-09 21:41:54 +02:00
|
|
|
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;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2019-08-07 20:16:50 +02:00
|
|
|
count = sizes_[0] * sizes_[1] * sizes_[2];
|
2019-05-09 21:41:54 +02:00
|
|
|
|
2019-08-07 20:16:50 +02:00
|
|
|
#ifdef USE_MPI
|
2019-11-14 15:36:39 +01:00
|
|
|
if( bdistributed ){
|
|
|
|
double globsum1{0.0};
|
|
|
|
size_t globcount{0};
|
2019-08-07 20:16:50 +02:00
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
MPI_Allreduce(reinterpret_cast<const void *>(&sum1),
|
|
|
|
reinterpret_cast<void *>(&globsum1),
|
|
|
|
1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
|
2019-08-07 20:16:50 +02:00
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
MPI_Allreduce(reinterpret_cast<const void *>(&count),
|
|
|
|
reinterpret_cast<void *>(&globcount),
|
|
|
|
1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, MPI_COMM_WORLD);
|
2019-08-07 20:16:50 +02:00
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
sum1 = globsum1;
|
|
|
|
count = globcount;
|
|
|
|
}
|
2019-10-09 17:25:15 +02:00
|
|
|
#endif
|
2019-08-07 20:16:50 +02:00
|
|
|
|
|
|
|
sum1 /= count;
|
2019-05-09 21:41:54 +02:00
|
|
|
|
|
|
|
return sum1;
|
|
|
|
}
|
|
|
|
|
|
|
|
template <typename functional, typename grid_t>
|
|
|
|
void assign_function_of_grids_r(const functional &f, const grid_t &g)
|
|
|
|
{
|
2019-11-14 15:36:39 +01:00
|
|
|
assert(g.size(0) == size(0) && g.size(1) == size(1));
|
2019-05-09 21:41:54 +02:00
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
#pragma omp parallel for
|
2019-05-09 21:41:54 +02:00
|
|
|
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);
|
|
|
|
|
|
|
|
elem = f(elemg);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
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)
|
|
|
|
{
|
2019-11-14 15:36:39 +01:00
|
|
|
assert(g1.size(0) == size(0) && g1.size(1) == size(1));
|
|
|
|
assert(g2.size(0) == size(0) && g2.size(1) == size(1));
|
2019-05-09 21:41:54 +02:00
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
#pragma omp parallel for
|
2019-05-09 21:41:54 +02:00
|
|
|
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);
|
|
|
|
|
|
|
|
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));
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
#pragma omp parallel for
|
2019-05-09 21:41:54 +02:00
|
|
|
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);
|
|
|
|
|
|
|
|
elem = f(elemg1, elemg2, elemg3);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2019-08-05 15:04:50 +02:00
|
|
|
template <typename functional, typename grid_t>
|
|
|
|
void assign_function_of_grids_k(const functional &f, const grid_t &g)
|
|
|
|
{
|
|
|
|
assert(g.size(0) == size(0) && g.size(1) == size(1)); // && g.size(2) == size(2) );
|
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
#pragma omp parallel for
|
2019-08-05 15:04:50 +02:00
|
|
|
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);
|
|
|
|
const auto &elemg = g.kelem(i, j, k);
|
|
|
|
|
|
|
|
elem = f(elemg);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2019-09-16 17:55:42 +02:00
|
|
|
template <typename functional, typename grid1_t, typename grid2_t>
|
|
|
|
void assign_function_of_grids_k(const functional &f, const grid1_t &g1, const grid2_t &g2)
|
|
|
|
{
|
|
|
|
assert(g1.size(0) == size(0) && g1.size(1) == size(1)); // && g.size(2) == size(2) );
|
|
|
|
assert(g2.size(0) == size(0) && g2.size(1) == size(1)); // && g.size(2) == size(2) );
|
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
#pragma omp parallel for
|
2019-09-16 17:55:42 +02:00
|
|
|
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);
|
|
|
|
const auto &elemg1 = g1.kelem(i, j, k);
|
|
|
|
const auto &elemg2 = g2.kelem(i, j, k);
|
|
|
|
|
2019-10-09 17:25:15 +02:00
|
|
|
elem = f(elemg1, elemg2);
|
2019-09-16 17:55:42 +02:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2019-11-03 15:54:17 +01:00
|
|
|
template <typename functional, typename grid_t>
|
|
|
|
void assign_function_of_grids_kdep(const functional &f, const grid_t &g)
|
|
|
|
{
|
|
|
|
assert(g.size(0) == size(0) && g.size(1) == size(1)); // && g.size(2) == size(2) );
|
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
#pragma omp parallel for
|
2019-11-03 15:54:17 +01:00
|
|
|
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);
|
|
|
|
const auto &elemg = g.kelem(i, j, k);
|
|
|
|
|
|
|
|
elem = f(this->get_k<real_t>(i, j, k), elemg);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2019-09-16 17:55:42 +02:00
|
|
|
template <typename functional, typename grid1_t, typename grid2_t>
|
|
|
|
void assign_function_of_grids_kdep(const functional &f, const grid1_t &g1, const grid2_t &g2)
|
|
|
|
{
|
2019-11-04 00:25:45 +01:00
|
|
|
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) );
|
2019-09-16 17:55:42 +02:00
|
|
|
|
2019-11-14 15:36:39 +01:00
|
|
|
#pragma omp parallel for
|
2019-11-04 00:25:45 +01:00
|
|
|
for (size_t i = 0; i < size(0); ++i)
|
2019-09-16 17:55:42 +02:00
|
|
|
{
|
2019-11-04 00:25:45 +01:00
|
|
|
for (size_t j = 0; j < size(1); ++j)
|
2019-09-16 17:55:42 +02:00
|
|
|
{
|
2019-11-04 00:25:45 +01:00
|
|
|
for (size_t k = 0; k < size(2); ++k)
|
2019-09-16 17:55:42 +02:00
|
|
|
{
|
|
|
|
auto &elem = this->kelem(i, j, k);
|
|
|
|
const auto &elemg1 = g1.kelem(i, j, k);
|
|
|
|
const auto &elemg2 = g2.kelem(i, j, k);
|
|
|
|
|
2019-10-09 17:25:15 +02:00
|
|
|
elem = f(this->get_k<real_t>(i, j, k), elemg1, elemg2);
|
2019-09-16 17:55:42 +02:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
template <typename functional>
|
|
|
|
void apply_function_k_dep(const functional &f)
|
|
|
|
{
|
2019-11-14 15:36:39 +01:00
|
|
|
#pragma omp parallel for
|
2019-05-09 21:41:54 +02:00
|
|
|
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>
|
|
|
|
void apply_function_r_dep(const functional &f)
|
|
|
|
{
|
2019-11-14 15:36:39 +01:00
|
|
|
#pragma omp parallel for
|
2019-05-09 21:41:54 +02:00
|
|
|
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 FourierTransformBackward(bool do_transform = true);
|
|
|
|
|
|
|
|
void FourierTransformForward(bool do_transform = true);
|
|
|
|
|
|
|
|
void ApplyNorm(void);
|
|
|
|
|
|
|
|
void FillRandomReal(unsigned long int seed = 123456ul);
|
|
|
|
|
2019-08-02 19:07:45 +02:00
|
|
|
void Write_to_HDF5(std::string fname, std::string datasetname) const;
|
2019-05-09 21:41:54 +02:00
|
|
|
|
2019-10-09 17:25:15 +02:00
|
|
|
void Write_PowerSpectrum(std::string ofname);
|
2019-05-09 21:41:54 +02:00
|
|
|
|
2019-10-09 17:25:15 +02:00
|
|
|
void Compute_PowerSpectrum(std::vector<double> &bin_k, std::vector<double> &bin_P, std::vector<double> &bin_eP, std::vector<size_t> &bin_count);
|
2019-05-19 12:05:04 +02:00
|
|
|
|
|
|
|
void Write_PDF(std::string ofname, int nbins = 1000, double scale = 1.0, double rhomin = 1e-3, double rhomax = 1e3);
|
2019-05-09 21:41:54 +02:00
|
|
|
|
2019-11-04 00:25:45 +01:00
|
|
|
void shift_field( const vec3<real_t>& s, bool transform_back=true )
|
2019-10-09 17:25:15 +02:00
|
|
|
{
|
2019-08-02 19:17:40 +02:00
|
|
|
FourierTransformForward();
|
|
|
|
apply_function_k_dep([&](auto x, auto k) -> ccomplex_t {
|
2019-11-14 15:36:39 +01:00
|
|
|
real_t shift;
|
|
|
|
if( bdistributed ){
|
|
|
|
shift = s.y * k[0] * get_dx()[0] + s.x * k[1] * get_dx()[1] + s.z * k[2] * get_dx()[2];
|
|
|
|
}else{
|
|
|
|
shift = s.x * k[0] * get_dx()[0] + s.y * k[1] * get_dx()[1] + s.z * k[2] * get_dx()[2];
|
|
|
|
}
|
2019-10-17 17:39:26 +02:00
|
|
|
return x * std::exp(ccomplex_t(0.0, shift));
|
2019-08-02 19:17:40 +02:00
|
|
|
});
|
2019-11-04 00:25:45 +01:00
|
|
|
if( transform_back ){
|
|
|
|
FourierTransformBackward();
|
|
|
|
}
|
2019-08-02 19:17:40 +02:00
|
|
|
}
|
|
|
|
|
2019-05-09 21:41:54 +02:00
|
|
|
void zero_DC_mode(void)
|
|
|
|
{
|
2019-10-09 17:25:15 +02:00
|
|
|
if (space_ == kspace_id)
|
|
|
|
{
|
2019-11-14 15:36:39 +01:00
|
|
|
if (CONFIG::MPI_task_rank == 0 || !bdistributed )
|
2019-10-09 17:25:15 +02:00
|
|
|
cdata_[0] = (data_t)0.0;
|
|
|
|
}
|
|
|
|
else
|
|
|
|
{
|
2019-05-12 17:39:15 +02:00
|
|
|
data_t sum = 0.0;
|
2019-08-02 19:17:40 +02:00
|
|
|
// #pragma omp parallel for reduction(+:sum)
|
2019-05-12 17:39:15 +02:00
|
|
|
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)
|
|
|
|
{
|
|
|
|
sum += this->relem(i, j, k);
|
2019-05-12 18:28:37 +02:00
|
|
|
}
|
2019-05-12 17:39:15 +02:00
|
|
|
}
|
|
|
|
}
|
2019-11-14 15:36:39 +01:00
|
|
|
if( bdistributed ){
|
2019-10-09 17:25:15 +02:00
|
|
|
#if defined(USE_MPI)
|
2019-11-14 15:36:39 +01:00
|
|
|
data_t glob_sum = 0.0;
|
|
|
|
MPI_Allreduce(reinterpret_cast<void *>(&sum), reinterpret_cast<void *>(&glob_sum),
|
|
|
|
1, GetMPIDatatype<data_t>(), MPI_SUM, MPI_COMM_WORLD);
|
|
|
|
sum = glob_sum;
|
2019-10-09 17:25:15 +02:00
|
|
|
#endif
|
2019-11-14 15:36:39 +01:00
|
|
|
}
|
2019-10-09 17:25:15 +02:00
|
|
|
sum /= sizes_[0] * sizes_[1] * sizes_[2];
|
2019-05-12 17:39:15 +02:00
|
|
|
|
2019-08-05 15:04:50 +02:00
|
|
|
#pragma omp parallel for
|
2019-05-12 17:39:15 +02:00
|
|
|
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)
|
|
|
|
{
|
|
|
|
this->relem(i, j, k) -= sum;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2019-05-09 21:41:54 +02:00
|
|
|
}
|
2019-08-05 15:04:50 +02:00
|
|
|
|
2019-10-09 17:25:15 +02:00
|
|
|
void dealias(void)
|
|
|
|
{
|
|
|
|
static const real_t kmax[3] =
|
|
|
|
{(real_t)(2.0 / 3.0 * (real_t)n_[0] / 2 * kfac_[0]),
|
|
|
|
(real_t)(2.0 / 3.0 * (real_t)n_[1] / 2 * kfac_[1]),
|
|
|
|
(real_t)(2.0 / 3.0 * (real_t)n_[2] / 2 * kfac_[2])};
|
|
|
|
|
|
|
|
//static const real_t kmax2 = kmax*kmax;
|
|
|
|
|
|
|
|
for (size_t i = 0; i < this->size(0); ++i)
|
|
|
|
{
|
|
|
|
for (size_t j = 0; j < this->size(1); ++j)
|
|
|
|
{
|
|
|
|
// size_t idx = (i * this->size(1) + j) * this->size(3);
|
|
|
|
for (size_t k = 0; k < this->size(2); ++k)
|
|
|
|
{
|
|
|
|
auto kk = get_k<real_t>(i, j, k);
|
|
|
|
//if (std::abs(kk[0]) > kmax[0] || std::abs(kk[1]) > kmax[1] || std::abs(kk[2]) > kmax[2])
|
|
|
|
if( kk.norm() > kmax[0] )
|
|
|
|
this->kelem(i,j,k) = 0.0;
|
|
|
|
// ++idx;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2019-05-07 01:05:16 +02:00
|
|
|
};
|