mirror of
https://github.com/cosmo-sims/monofonIC.git
synced 2024-09-19 17:03:45 +02:00
cleaned up PANPHASIA plugin so that now there is only one PANPHASIA and it selects automatically if it is 'ho' or v1
This commit is contained in:
Oliver Hahn
parent
e3f5d65e0a
commit
d37cfc79b6
5 changed files with 118 additions and 672 deletions
|
|
@ -150,7 +150,7 @@ endif()
|
|||
endif(ENABLE_PANPHASIA)
|
||||
########################################################################################################################
|
||||
# PANPHASIA HO (High-Order, new version)
|
||||
option(ENABLE_PANPHASIA_HO "Enable PANPHASIA-HO random number generator" ON)
|
||||
# option(ENABLE_PANPHASIA_HO "Enable PANPHASIA-HO random number generator" ON)
|
||||
########################################################################################################################
|
||||
# INCLUDES
|
||||
include_directories(${PROJECT_SOURCE_DIR}/include ${PROJECT_SOURCE_DIR}/external)
|
||||
|
|
@ -173,12 +173,7 @@ if(ENABLE_PANPHASIA)
|
|||
list (APPEND SOURCES
|
||||
${PROJECT_SOURCE_DIR}/external/panphasia/panphasia_routines.f
|
||||
${PROJECT_SOURCE_DIR}/external/panphasia/generic_lecuyer.f90
|
||||
)
|
||||
endif()
|
||||
|
||||
if(ENABLE_PANPHASIA_HO)
|
||||
list (APPEND SOURCES
|
||||
${PROJECT_SOURCE_DIR}/external/panphasia_ho/main.c
|
||||
#${PROJECT_SOURCE_DIR}/external/panphasia_ho/main.c
|
||||
${PROJECT_SOURCE_DIR}/external/panphasia_ho/high_order_panphasia_routines.c
|
||||
${PROJECT_SOURCE_DIR}/external/panphasia_ho/pan_mpi_routines.c
|
||||
${PROJECT_SOURCE_DIR}/external/panphasia_ho/uniform_rand_threefry4x64.c
|
||||
|
|
@ -259,9 +254,9 @@ if(ENABLE_PANPHASIA)
|
|||
target_compile_definitions(${PRGNAME} PRIVATE "USE_PANPHASIA")
|
||||
endif(ENABLE_PANPHASIA)
|
||||
|
||||
if(ENABLE_PANPHASIA_HO)
|
||||
target_compile_definitions(${PRGNAME} PRIVATE "USE_PANPHASIA_HO")
|
||||
endif(ENABLE_PANPHASIA_HO)
|
||||
# if(ENABLE_PANPHASIA_HO)
|
||||
# target_compile_definitions(${PRGNAME} PRIVATE "USE_PANPHASIA_HO")
|
||||
# endif(ENABLE_PANPHASIA_HO)
|
||||
|
||||
if(ENABLE_PLT)
|
||||
target_compile_definitions(${PRGNAME} PRIVATE "ENABLE_PLT")
|
||||
|
|
|
|||
|
|
@ -16,7 +16,7 @@ LPTorder = 3 # order of the LPT to be used (1,2 or 3)
|
|||
DoBaryons = no # also do baryon ICs?
|
||||
DoBaryonVrel = no # if doing baryons, incl. also relative velocity to linear order?
|
||||
|
||||
DoFixing = yes # do mode fixing à la Angulo&Pontzen (https://arxiv.org/abs/1603.05253)
|
||||
DoFixing = no # do mode fixing à la Angulo&Pontzen (https://arxiv.org/abs/1603.05253)
|
||||
DoInversion = no # invert phases (for paired simulations)
|
||||
|
||||
ParticleLoad = sc # particle load, can be 'sc' (1x), 'bcc' (2x) or 'fcc' (4x)
|
||||
|
|
|
|||
|
|
@ -136,16 +136,14 @@ inline void HDFGetDatasetExtent( const std::string Filename, const std::string O
|
|||
|
||||
int ndims = H5Sget_simple_extent_ndims( HDF_DataspaceID );
|
||||
|
||||
hsize_t *dimsize = new hsize_t[ndims];
|
||||
std::vector<hsize_t> dimsize(ndims,0);
|
||||
|
||||
H5Sget_simple_extent_dims( HDF_DataspaceID, dimsize, NULL );
|
||||
H5Sget_simple_extent_dims( HDF_DataspaceID, &dimsize[0], NULL );
|
||||
|
||||
Extent.clear();
|
||||
for(int i=0; i<ndims; ++i )
|
||||
Extent.push_back( dimsize[i] );
|
||||
|
||||
delete[] dimsize;
|
||||
|
||||
H5Sclose( HDF_DataspaceID );
|
||||
H5Dclose( HDF_DatasetID );
|
||||
H5Fclose( HDF_FileID );
|
||||
|
|
@ -323,8 +321,8 @@ inline void HDFReadVectorSelect( const std::string Filename, const std::string O
|
|||
//... get space associated with dataset and its extensions
|
||||
HDF_DataspaceID = H5Dget_space( HDF_DatasetID );
|
||||
int ndims = H5Sget_simple_extent_ndims( HDF_DataspaceID );
|
||||
hsize_t dimsize[ndims];
|
||||
H5Sget_simple_extent_dims( HDF_DataspaceID, dimsize, NULL );
|
||||
std::vector<hsize_t> dimsize(ndims,0);
|
||||
H5Sget_simple_extent_dims( HDF_DataspaceID, &dimsize[0], NULL );
|
||||
|
||||
hsize_t block[2];
|
||||
block[0] = ii.size();
|
||||
|
|
@ -601,9 +599,9 @@ inline void HDFReadGroupAttribute( const std::string Filename, const std::string
|
|||
|
||||
int ndims = H5Sget_simple_extent_ndims( HDF_DataspaceID );
|
||||
|
||||
hsize_t dimsize[ndims];
|
||||
std::vector<hsize_t> dimsize(ndims,0);
|
||||
|
||||
H5Sget_simple_extent_dims( HDF_DataspaceID, dimsize, NULL );
|
||||
H5Sget_simple_extent_dims( HDF_DataspaceID, &dimsize[0], NULL );
|
||||
|
||||
HDF_StorageSize = 1;
|
||||
for(int i=0; i<ndims; ++i )
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
// This file is part of monofonIC (MUSIC2)
|
||||
// A software package to generate ICs for cosmological simulations
|
||||
// Copyright (C) 2020 by Oliver Hahn and Adrian Jenkins (this file)
|
||||
// Copyright (C) 2021 by Oliver Hahn and Adrian Jenkins (this file)
|
||||
// but see distinct licensing for PANPHASIA below
|
||||
//
|
||||
// monofonIC is free software: you can redistribute it and/or modify
|
||||
|
|
@ -39,108 +39,16 @@
|
|||
|
||||
#include <grid_fft.hh>
|
||||
|
||||
const int maxdim = 60, maxlev = 50, maxpow = 3 * maxdim;
|
||||
typedef int rand_offset_[5];
|
||||
typedef struct
|
||||
namespace PANPHASIA2
|
||||
{
|
||||
int state[133]; // Nstore = Nstate (=5) + Nbatch (=128)
|
||||
int need_fill;
|
||||
int pos;
|
||||
} rand_state_;
|
||||
|
||||
/* pan_state_ struct -- corresponds to respective fortran module in panphasia_routines.f
|
||||
* data structure that contains all panphasia state variables
|
||||
* it needs to get passed between the fortran routines to enable
|
||||
* thread-safe execution.
|
||||
*/
|
||||
typedef struct
|
||||
{
|
||||
int base_state[5], base_lev_start[5][maxdim + 1];
|
||||
rand_offset_ poweroffset[maxpow + 1], superjump;
|
||||
rand_state_ current_state[maxpow + 2];
|
||||
|
||||
int layer_min, layer_max, indep_field;
|
||||
|
||||
long long xorigin_store[2][2][2], yorigin_store[2][2][2], zorigin_store[2][2][2];
|
||||
int lev_common, layer_min_store, layer_max_store;
|
||||
long long ix_abs_store, iy_abs_store, iz_abs_store, ix_per_store, iy_per_store, iz_per_store, ix_rel_store,
|
||||
iy_rel_store, iz_rel_store;
|
||||
double exp_coeffs[8][8][maxdim + 2];
|
||||
long long xcursor[maxdim + 1], ycursor[maxdim + 1], zcursor[maxdim + 1];
|
||||
int ixshift[2][2][2], iyshift[2][2][2], izshift[2][2][2];
|
||||
|
||||
double cell_data[9][8];
|
||||
int ixh_last, iyh_last, izh_last;
|
||||
int init;
|
||||
|
||||
int init_cell_props;
|
||||
int init_lecuyer_state;
|
||||
long long p_xcursor[62], p_ycursor[62], p_zcursor[62];
|
||||
|
||||
} pan_state_;
|
||||
|
||||
extern "C"
|
||||
{
|
||||
void start_panphasia_(pan_state_ *lstate, const char *descriptor, int *ngrid, int *bverbose);
|
||||
|
||||
void parse_descriptor_(const char *descriptor, int16_t *l, int32_t *ix, int32_t *iy, int32_t *iz, int16_t *side1,
|
||||
int16_t *side2, int16_t *side3, int32_t *check_int, char *name);
|
||||
|
||||
void panphasia_cell_properties_(pan_state_ *lstate, int *ixcell, int *iycell, int *izcell, double *cell_prop);
|
||||
|
||||
void adv_panphasia_cell_properties_(pan_state_ *lstate, int *ixcell, int *iycell, int *izcell, int *layer_min,
|
||||
int *layer_max, int *indep_field, double *cell_prop);
|
||||
|
||||
void set_phases_and_rel_origin_(pan_state_ *lstate, const char *descriptor, int *lev, long long *ix_rel,
|
||||
long long *iy_rel, long long *iz_rel, int *VERBOSE);
|
||||
}
|
||||
|
||||
struct panphasia_descriptor
|
||||
{
|
||||
int16_t wn_level_base;
|
||||
int32_t i_xorigin_base, i_yorigin_base, i_zorigin_base;
|
||||
int16_t i_base, i_base_y, i_base_z;
|
||||
int32_t check_rand;
|
||||
std::string name;
|
||||
|
||||
explicit panphasia_descriptor(std::string dstring)
|
||||
extern "C"
|
||||
{
|
||||
char tmp[100];
|
||||
std::memset(tmp, ' ', 100);
|
||||
parse_descriptor_(dstring.c_str(), &wn_level_base, &i_xorigin_base, &i_yorigin_base, &i_zorigin_base, &i_base,
|
||||
&i_base_y, &i_base_z, &check_rand, tmp);
|
||||
for (int i = 0; i < 100; i++)
|
||||
if (tmp[i] == ' ')
|
||||
{
|
||||
tmp[i] = '\0';
|
||||
break;
|
||||
}
|
||||
name = tmp;
|
||||
name.erase(std::remove(name.begin(), name.end(), ' '), name.end());
|
||||
#include <panphasia_ho/panphasia_functions.h>
|
||||
extern size_t descriptor_base_size;
|
||||
}
|
||||
};
|
||||
|
||||
// greatest common divisor
|
||||
int gcd(int a, int b)
|
||||
{
|
||||
if (b == 0)
|
||||
return a;
|
||||
return gcd(b, a % b);
|
||||
}
|
||||
|
||||
// least common multiple
|
||||
int lcm(int a, int b) { return abs(a * b) / gcd(a, b); }
|
||||
|
||||
// Two or largest power of 2 less than the argument
|
||||
int largest_power_two_lte(int b)
|
||||
{
|
||||
int a = 1;
|
||||
if (b <= a)
|
||||
return a;
|
||||
while (2 * a < b)
|
||||
a = 2 * a;
|
||||
return a;
|
||||
}
|
||||
#include "PANPHASIA.hh"
|
||||
|
||||
class RNG_panphasia : public RNG_plugin
|
||||
{
|
||||
|
|
@ -148,63 +56,15 @@ private:
|
|||
protected:
|
||||
std::string descriptor_string_;
|
||||
int num_threads_;
|
||||
int levelmin_, levelmin_final_, levelmax_, ngrid_, ngrid_panphasia_;
|
||||
bool incongruent_fields_;
|
||||
double inter_grid_phase_adjustment_;
|
||||
// double translation_phase_;
|
||||
pan_state_ *lstate;
|
||||
int grid_p_;
|
||||
int coordinate_system_shift_[3];
|
||||
int ix_abs_[3], ix_per_[3], ix_rel_[3], level_p_, lextra_;
|
||||
int panphasia_mode_;
|
||||
size_t grid_res_;
|
||||
real_t boxlength_;
|
||||
|
||||
void clear_panphasia_thread_states(void)
|
||||
{
|
||||
for (int i = 0; i < num_threads_; ++i)
|
||||
{
|
||||
lstate[i].init = 0;
|
||||
lstate[i].init_cell_props = 0;
|
||||
lstate[i].init_lecuyer_state = 0;
|
||||
}
|
||||
}
|
||||
|
||||
void initialize_for_grid_structure(void)
|
||||
{
|
||||
// if ngrid is not a multiple of i_base, then we need to enlarge and then sample down
|
||||
ngrid_ = pcf_->get_value<size_t>("setup", "GridRes");
|
||||
int ngridminsize_panphasia = pcf_->get_value_safe<int>("random", "PanphasiaMinRootResolution",512);
|
||||
|
||||
grid_p_ = pdescriptor_->i_base;
|
||||
|
||||
lextra_ = (log10((double)ngrid_ / (double)grid_p_) + 0.001) / log10(2.0);
|
||||
// lmin
|
||||
|
||||
ngrid_panphasia_ = (1 << lextra_) * grid_p_;
|
||||
|
||||
while( ngrid_panphasia_ < ngridminsize_panphasia ){
|
||||
lextra_++;
|
||||
ngrid_panphasia_*=2;
|
||||
}
|
||||
assert( ngrid_panphasia_ >= ngridminsize_panphasia);
|
||||
|
||||
|
||||
clear_panphasia_thread_states();
|
||||
|
||||
music::ilog.Print("PANPHASIA: using grid size %lld (level=%d)",ngrid_panphasia_, lextra_);
|
||||
if (ngridminsize_panphasia<512)
|
||||
music::ilog.Print("PANPHASIA WARNING: PanphasiaMinRootResolution = %d below minimum recommended of 512",ngridminsize_panphasia);
|
||||
music::ilog.Print("PANPHASIA: running with %d threads", num_threads_, ngrid_panphasia_ );
|
||||
|
||||
coordinate_system_shift_[0] = -pcf_->get_value_safe<int>("setup", "shift_x", 0);
|
||||
coordinate_system_shift_[1] = -pcf_->get_value_safe<int>("setup", "shift_y", 0);
|
||||
coordinate_system_shift_[2] = -pcf_->get_value_safe<int>("setup", "shift_z", 0);
|
||||
}
|
||||
|
||||
std::unique_ptr<panphasia_descriptor> pdescriptor_;
|
||||
PANPHASIA1::RNG *ppan1_rng_;
|
||||
|
||||
public:
|
||||
explicit RNG_panphasia(config_file &cf) : RNG_plugin(cf)
|
||||
{
|
||||
descriptor_string_ = pcf_->get_value<std::string>("random", "descriptor");
|
||||
|
||||
#ifdef _OPENMP
|
||||
num_threads_ = omp_get_max_threads();
|
||||
|
|
@ -212,334 +72,116 @@ public:
|
|||
num_threads_ = 1;
|
||||
#endif
|
||||
|
||||
// create independent state descriptions for each thread
|
||||
lstate = new pan_state_[num_threads_];
|
||||
descriptor_string_ = pcf_->get_value<std::string>("random", "descriptor");
|
||||
grid_res_ = pcf_->get_value<size_t>("setup", "GridRes");
|
||||
boxlength_ = pcf_->get_value<real_t>("setup", "BoxLength");
|
||||
|
||||
// parse the descriptor for its properties
|
||||
pdescriptor_ = std::make_unique<panphasia_descriptor>(descriptor_string_);
|
||||
if( pcf_->get_value<bool>("setup", "DoFixing") ){
|
||||
music::flog << "Fixing all the modes to the mean power negates any advantage of using the Panphasia field.\n";
|
||||
music::flog << "With the new panphasia ho it is possible by choosing the descriptor to fix the largest modes without losing the ability to resimulate to much higher resolution.\n";
|
||||
throw std::runtime_error("PANPHASIA: incompatible parameter.");
|
||||
}
|
||||
|
||||
music::ilog.Print("PANPHASIA: descriptor \'%s\' is base %d,", pdescriptor_->name.c_str(), pdescriptor_->i_base);
|
||||
panphasia_mode_ = 0;
|
||||
PANPHASIA2::parse_and_validate_descriptor_(descriptor_string_.c_str(), &panphasia_mode_);
|
||||
|
||||
// write panphasia base size into config file for the grid construction
|
||||
// as the gridding unit we use the least common multiple of 2 and i_base
|
||||
std::stringstream ss;
|
||||
//ARJ ss << lcm(2, pdescriptor_->i_base);
|
||||
//ss << two_or_largest_power_two_less_than(pdescriptor_->i_base);//ARJ
|
||||
ss << 2; //ARJ - set gridding unit to two
|
||||
pcf_->insert_value("setup", "gridding_unit", ss.str());
|
||||
ss.str(std::string());
|
||||
ss << pdescriptor_->i_base;
|
||||
pcf_->insert_value("random", "base_unit", ss.str());
|
||||
|
||||
this->initialize_for_grid_structure();
|
||||
if (panphasia_mode_ == 0)
|
||||
{
|
||||
ppan1_rng_ = new PANPHASIA1::RNG(&cf);
|
||||
}
|
||||
}
|
||||
|
||||
~RNG_panphasia() { delete[] lstate; }
|
||||
~RNG_panphasia()
|
||||
{
|
||||
if (panphasia_mode_ == 0)
|
||||
{
|
||||
delete ppan1_rng_;
|
||||
}
|
||||
}
|
||||
|
||||
bool isMultiscale() const { return true; }
|
||||
|
||||
void Run_Panphasia_Highorder(Grid_FFT<real_t> &g);
|
||||
|
||||
void Fill_Grid(Grid_FFT<real_t> &g)
|
||||
{
|
||||
auto sinc = [](real_t x) { return (std::fabs(x) > 1e-16) ? std::sin(x) / x : 1.0; };
|
||||
auto dsinc = [](real_t x) { return (std::fabs(x) > 1e-16) ? (x * std::cos(x) - std::sin(x)) / (x * x) : 0.0; };
|
||||
const real_t sqrt3{std::sqrt(3.0)}, sqrt27{std::sqrt(27.0)};
|
||||
|
||||
// we will overwrite 'g', we can deallocate it while we prepare the panphasia field
|
||||
g.reset();
|
||||
|
||||
clear_panphasia_thread_states();
|
||||
|
||||
// temporaries
|
||||
Grid_FFT<real_t> g0({size_t(ngrid_panphasia_),size_t(ngrid_panphasia_),size_t(ngrid_panphasia_)}, g.length_);
|
||||
Grid_FFT<real_t> g1({size_t(ngrid_panphasia_),size_t(ngrid_panphasia_),size_t(ngrid_panphasia_)}, g.length_);
|
||||
Grid_FFT<real_t> g2({size_t(ngrid_panphasia_),size_t(ngrid_panphasia_),size_t(ngrid_panphasia_)}, g.length_);
|
||||
Grid_FFT<real_t> g3({size_t(ngrid_panphasia_),size_t(ngrid_panphasia_),size_t(ngrid_panphasia_)}, g.length_);
|
||||
Grid_FFT<real_t> g4({size_t(ngrid_panphasia_),size_t(ngrid_panphasia_),size_t(ngrid_panphasia_)}, g.length_);
|
||||
|
||||
double t1 = get_wtime();
|
||||
// double tp = t1;
|
||||
|
||||
#pragma omp parallel
|
||||
switch (panphasia_mode_)
|
||||
{
|
||||
#ifdef _OPENMP
|
||||
const int mythread = omp_get_thread_num();
|
||||
#else
|
||||
const int mythread = 0;
|
||||
#endif
|
||||
|
||||
//int odd_x, odd_y, odd_z;
|
||||
//int ng_level = ngrid_ * (1 << (level - levelmin_)); // full resolution of current level
|
||||
case 0: // old mode
|
||||
music::ilog << "PANPHASIA: Old descriptor" << std::endl;
|
||||
ppan1_rng_->Fill(g);
|
||||
break;
|
||||
|
||||
int verbosity = (mythread == 0);
|
||||
char descriptor[100];
|
||||
std::memset(descriptor, 0, 100);
|
||||
std::memcpy(descriptor, descriptor_string_.c_str(), descriptor_string_.size());
|
||||
case 1: // PANPHASIA HO descriptor
|
||||
music::ilog << "PANPHASIA: New descriptor" << std::endl;
|
||||
this->Run_Panphasia_Highorder(g);
|
||||
break;
|
||||
|
||||
start_panphasia_(&lstate[mythread], descriptor, &ngrid_panphasia_, &verbosity);
|
||||
|
||||
{
|
||||
panphasia_descriptor d(descriptor_string_);
|
||||
|
||||
int level_p = d.wn_level_base + lextra_;
|
||||
|
||||
lstate[mythread].layer_min = 0;
|
||||
lstate[mythread].layer_max = level_p;
|
||||
lstate[mythread].indep_field = 1;
|
||||
|
||||
long long ix_rel[3];
|
||||
ix_rel[0] = 0; //ileft_corner_p[0];
|
||||
ix_rel[1] = 0; //ileft_corner_p[1];
|
||||
ix_rel[2] = 0; //ileft_corner_p[2];
|
||||
|
||||
set_phases_and_rel_origin_(&lstate[mythread], descriptor, &level_p, &ix_rel[0], &ix_rel[1], &ix_rel[2],
|
||||
&verbosity);
|
||||
}
|
||||
|
||||
if (verbosity)
|
||||
t1 = get_wtime();
|
||||
|
||||
std::array<double, 9> cell_prop;
|
||||
pan_state_ *ps = &lstate[mythread];
|
||||
|
||||
#pragma omp for //nowait
|
||||
for (size_t i = 0; i < g0.size(0); i += 2)
|
||||
{
|
||||
const int ixmax(std::min<int>(2,g0.size(0)-i));
|
||||
for (size_t j = 0; j < g0.size(1); j += 2)
|
||||
{
|
||||
const int iymax(std::min<int>(2,g0.size(1)-j));
|
||||
for (size_t k = 0; k < g0.size(2); k += 2)
|
||||
{
|
||||
const int izmax(std::min<int>(2,g0.size(2)-k));
|
||||
|
||||
// ARJ - added inner set of loops to speed up evaluation of Panphasia
|
||||
for (int ix = 0; ix < ixmax; ++ix)
|
||||
{
|
||||
for (int iy = 0; iy < iymax; ++iy)
|
||||
{
|
||||
for (int iz = 0; iz < izmax; ++iz)
|
||||
{
|
||||
int ilocal = i + ix;
|
||||
int jlocal = j + iy;
|
||||
int klocal = k + iz;
|
||||
|
||||
int iglobal = ilocal + g0.local_0_start_;
|
||||
int jglobal = jlocal;
|
||||
int kglobal = klocal;
|
||||
|
||||
adv_panphasia_cell_properties_(ps, &iglobal, &jglobal, &kglobal, &ps->layer_min,
|
||||
&ps->layer_max, &ps->indep_field, &cell_prop[0]);
|
||||
|
||||
g0.relem(ilocal, jlocal, klocal) = cell_prop[0];
|
||||
g1.relem(ilocal, jlocal, klocal) = cell_prop[4];
|
||||
g2.relem(ilocal, jlocal, klocal) = cell_prop[2];
|
||||
g3.relem(ilocal, jlocal, klocal) = cell_prop[1];
|
||||
g4.relem(ilocal, jlocal, klocal) = cell_prop[8];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
} // end omp parallel region
|
||||
|
||||
g0.FourierTransformForward();
|
||||
g1.FourierTransformForward();
|
||||
g2.FourierTransformForward();
|
||||
g3.FourierTransformForward();
|
||||
g4.FourierTransformForward();
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < g0.size(0); i++)
|
||||
{
|
||||
for (size_t j = 0; j < g0.size(1); j++)
|
||||
{
|
||||
for (size_t k = 0; k < g0.size(2); k++)
|
||||
{
|
||||
if (!g0.is_nyquist_mode(i, j, k))
|
||||
{
|
||||
auto kvec = g0.get_k<real_t>(i, j, k);
|
||||
|
||||
auto argx = 0.5 * M_PI * kvec[0] / g0.kny_[0];
|
||||
auto argy = 0.5 * M_PI * kvec[1] / g0.kny_[1];
|
||||
auto argz = 0.5 * M_PI * kvec[2] / g0.kny_[2];
|
||||
|
||||
auto fx = real_t(sinc(argx));
|
||||
auto gx = ccomplex_t(0.0, dsinc(argx));
|
||||
auto fy = real_t(sinc(argy));
|
||||
auto gy = ccomplex_t(0.0, dsinc(argy));
|
||||
auto fz = real_t(sinc(argz));
|
||||
auto gz = ccomplex_t(0.0, dsinc(argz));
|
||||
|
||||
auto temp = (fx + sqrt3 * gx) * (fy + sqrt3 * gy) * (fz + sqrt3 * gz);
|
||||
auto magnitude = real_t(std::sqrt(1.0 - std::fabs(temp * temp)));
|
||||
|
||||
auto y0(g0.kelem(i, j, k)), y1(g1.kelem(i, j, k)), y2(g2.kelem(i, j, k)), y3(g3.kelem(i, j, k)), y4(g4.kelem(i, j, k));
|
||||
|
||||
g0.kelem(i, j, k) = y0 * fx * fy * fz
|
||||
+ sqrt3 * (y1 * gx * fy * fz + y2 * fx * gy * fz + y3 * fx * fy * gz)
|
||||
+ y4 * magnitude;
|
||||
}
|
||||
else
|
||||
{
|
||||
g0.kelem(i, j, k) = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
default: // unknown PANPHASIA mode
|
||||
music::elog << "PANPHASIA: Something went wrong with descriptor" << std::endl;
|
||||
abort();
|
||||
break;
|
||||
}
|
||||
|
||||
// music::ilog.Print("\033[31mtiming [build panphasia field]: %f s\033[0m", get_wtime() - tp);
|
||||
// tp = get_wtime();
|
||||
|
||||
g1.FourierTransformBackward(false);
|
||||
g2.FourierTransformBackward(false);
|
||||
g3.FourierTransformBackward(false);
|
||||
g4.FourierTransformBackward(false);
|
||||
|
||||
#pragma omp parallel
|
||||
{
|
||||
#ifdef _OPENMP
|
||||
const int mythread = omp_get_thread_num();
|
||||
#else
|
||||
const int mythread = 0;
|
||||
#endif
|
||||
|
||||
// int odd_x, odd_y, odd_z;
|
||||
int verbosity = (mythread == 0);
|
||||
char descriptor[100];
|
||||
std::memset(descriptor, 0, 100);
|
||||
std::memcpy(descriptor, descriptor_string_.c_str(), descriptor_string_.size());
|
||||
|
||||
start_panphasia_(&lstate[mythread], descriptor, &ngrid_panphasia_, &verbosity);
|
||||
|
||||
{
|
||||
panphasia_descriptor d(descriptor_string_);
|
||||
|
||||
int level_p = d.wn_level_base + lextra_;
|
||||
|
||||
lstate[mythread].layer_min = 0;
|
||||
lstate[mythread].layer_max = level_p;
|
||||
lstate[mythread].indep_field = 1;
|
||||
|
||||
long long ix_rel[3];
|
||||
ix_rel[0] = 0; //ileft_corner_p[0];
|
||||
ix_rel[1] = 0; //ileft_corner_p[1];
|
||||
ix_rel[2] = 0; //ileft_corner_p[2];
|
||||
|
||||
set_phases_and_rel_origin_(&lstate[mythread], descriptor, &level_p, &ix_rel[0], &ix_rel[1], &ix_rel[2],
|
||||
&verbosity);
|
||||
}
|
||||
|
||||
if (verbosity)
|
||||
t1 = get_wtime();
|
||||
|
||||
//***************************************************************
|
||||
// Process Panphasia values: p110, p011, p101, p111
|
||||
//****************************************************************
|
||||
std::array<double,9> cell_prop;
|
||||
pan_state_ *ps = &lstate[mythread];
|
||||
|
||||
#pragma omp for //nowait
|
||||
for (size_t i = 0; i < g1.size(0); i += 2)
|
||||
{
|
||||
const int ixmax(std::min<int>(2,g1.size(0)-i));
|
||||
for (size_t j = 0; j < g1.size(1); j += 2)
|
||||
{
|
||||
const int iymax(std::min<int>(2,g1.size(1)-j));
|
||||
for (size_t k = 0; k < g1.size(2); k += 2)
|
||||
{
|
||||
const int izmax(std::min<int>(2,g1.size(2)-k));
|
||||
|
||||
// ARJ - added inner set of loops to speed up evaluation of Panphasia
|
||||
for (int ix = 0; ix < ixmax; ++ix)
|
||||
{
|
||||
for (int iy = 0; iy < iymax; ++iy)
|
||||
{
|
||||
for (int iz = 0; iz < izmax; ++iz)
|
||||
{
|
||||
int ilocal = i + ix;
|
||||
int jlocal = j + iy;
|
||||
int klocal = k + iz;
|
||||
|
||||
int iglobal = ilocal + g1.local_0_start_;
|
||||
int jglobal = jlocal;
|
||||
int kglobal = klocal;
|
||||
|
||||
adv_panphasia_cell_properties_(ps, &iglobal, &jglobal, &kglobal, &ps->layer_min,
|
||||
&ps->layer_max, &ps->indep_field, &cell_prop[0]);
|
||||
|
||||
g1.relem(ilocal, jlocal, klocal) = cell_prop[6];
|
||||
g2.relem(ilocal, jlocal, klocal) = cell_prop[3];
|
||||
g3.relem(ilocal, jlocal, klocal) = cell_prop[5];
|
||||
g4.relem(ilocal, jlocal, klocal) = cell_prop[7];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
} // end omp parallel region
|
||||
|
||||
// music::ilog.Print("\033[31mtiming [adv_panphasia_cell_properties2]: %f s \033[0m", get_wtime() - tp);
|
||||
// tp = get_wtime();
|
||||
|
||||
/////////////////////////////////////////////////////////////////////////
|
||||
// transform and convolve with Legendres
|
||||
g1.FourierTransformForward();
|
||||
g2.FourierTransformForward();
|
||||
g3.FourierTransformForward();
|
||||
g4.FourierTransformForward();
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < g0.size(0); i++)
|
||||
{
|
||||
for (size_t j = 0; j < g0.size(1); j++)
|
||||
{
|
||||
for (size_t k = 0; k < g0.size(2); k++)
|
||||
{
|
||||
if (!g0.is_nyquist_mode(i, j, k))
|
||||
{
|
||||
auto kvec = g0.get_k<real_t>(i, j, k);
|
||||
|
||||
auto argx = 0.5 * M_PI * kvec[0] / g0.kny_[0];
|
||||
auto argy = 0.5 * M_PI * kvec[1] / g0.kny_[1];
|
||||
auto argz = 0.5 * M_PI * kvec[2] / g0.kny_[2];
|
||||
|
||||
auto fx = real_t(sinc(argx));
|
||||
auto gx = ccomplex_t(0.0, dsinc(argx));
|
||||
auto fy = real_t(sinc(argy));
|
||||
auto gy = ccomplex_t(0.0, dsinc(argy));
|
||||
auto fz = real_t(sinc(argz));
|
||||
auto gz = ccomplex_t(0.0, dsinc(argz));
|
||||
|
||||
auto y1(g1.kelem(i, j, k)), y2(g2.kelem(i, j, k)), y3(g3.kelem(i, j, k)), y4(g4.kelem(i, j, k));
|
||||
|
||||
g0.kelem(i, j, k) += real_t(3.0) * (y1 * gx * gy * fz + y2 * fx * gy * gz + y3 * gx * fy * gz) + sqrt27 * y4 * gx * gy * gz;
|
||||
|
||||
// do final phase shift to account for corner centered coordinates vs. cell centers
|
||||
auto phase_shift = - 0.5 * M_PI * ( kvec[0] / g0.kny_[0]
|
||||
+ kvec[1] /g0.kny_[1] + kvec[2] / g0.kny_[2]);
|
||||
|
||||
g0.kelem(i, j, k) *= std::exp( ccomplex_t(0,phase_shift) );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
g1.reset();
|
||||
g2.reset();
|
||||
g3.reset();
|
||||
g4.reset();
|
||||
|
||||
g.allocate();
|
||||
g0.FourierInterpolateCopyTo( g );
|
||||
|
||||
music::ilog.Print("time for calculating PANPHASIA field : %f s, %f µs/cell", get_wtime() - t1,
|
||||
1e6 * (get_wtime() - t1) / g.global_size(0) / g.global_size(1) / g.global_size(2));
|
||||
music::ilog.Print("PANPHASIA k-space statistices: mean Re = %f, std = %f", g.mean(), g.std());
|
||||
}
|
||||
};
|
||||
|
||||
void RNG_panphasia::Run_Panphasia_Highorder(Grid_FFT<real_t> &g)
|
||||
{
|
||||
int verbose = 0;
|
||||
int error;
|
||||
size_t x0 = 0, y0 = 0, z0 = 0;
|
||||
size_t rel_level;
|
||||
int fdim = 2; //Option to scale Fourier grid dimension relative to Panphasia coefficient grid
|
||||
|
||||
//char descriptor[300] = "[Panph6,L20,(424060,82570,148256),S1,KK0,CH-999,Auriga_100_vol2]";
|
||||
|
||||
PANPHASIA2::PANPHASIA_init_descriptor_(descriptor_string_.c_str(), &verbose);
|
||||
|
||||
printf("Descriptor %s\n ngrid_load %lu\n", descriptor_string_.c_str(), grid_res_);
|
||||
|
||||
// Choose smallest value of level to equal of exceed grid_res_)
|
||||
|
||||
for (rel_level = 0; fdim * (PANPHASIA2::descriptor_base_size <<rel_level) < grid_res_; rel_level++);
|
||||
|
||||
printf("Setting relative level = %lu\n", rel_level);
|
||||
|
||||
if ((error = PANPHASIA2::PANPHASIA_init_level_(&rel_level, &x0, &y0, &z0, &verbose)))
|
||||
{
|
||||
printf("Abort: PANPHASIA_init_level_ :error code %d\n", error);
|
||||
abort();
|
||||
};
|
||||
|
||||
//======================= FFTW ==============================
|
||||
|
||||
ptrdiff_t alloc_local, local_n0, local_0_start;
|
||||
|
||||
size_t N0 = fdim * (PANPHASIA2::descriptor_base_size << rel_level);
|
||||
|
||||
alloc_local = FFTW_MPI_LOCAL_SIZE_3D(N0, N0, N0 + 2, MPI_COMM_WORLD, &local_n0, &local_0_start);
|
||||
|
||||
Grid_FFT<real_t> pan_grid({{N0, N0, N0}}, {{boxlength_, boxlength_, boxlength_}});
|
||||
|
||||
assert(pan_grid.n_[0] == N0);
|
||||
assert(pan_grid.n_[1] == N0);
|
||||
assert(pan_grid.n_[2] == N0);
|
||||
assert(pan_grid.local_0_start_ == local_0_start);
|
||||
assert(pan_grid.local_0_size_ == local_n0);
|
||||
assert(pan_grid.ntot_ == size_t(alloc_local));
|
||||
|
||||
pan_grid.FourierTransformForward(false);
|
||||
|
||||
FFTW_COMPLEX *Panphasia_White_Noise_Field = reinterpret_cast<FFTW_COMPLEX *>(&pan_grid.data_[0]);
|
||||
|
||||
// Panphasia_White_Noise_Field = FFTW_ALLOC_COMPLEX(alloc_local);
|
||||
|
||||
if ((error = PANPHASIA2::PANPHASIA_compute_kspace_field_(rel_level, N0, local_n0, local_0_start, Panphasia_White_Noise_Field)))
|
||||
{
|
||||
music::elog << "Error code from PANPHASIA_compute ... (ErrCode = " << error << ")" << std::endl;
|
||||
};
|
||||
|
||||
pan_grid.FourierInterpolateCopyTo(g);
|
||||
}
|
||||
namespace
|
||||
{
|
||||
RNG_plugin_creator_concrete<RNG_panphasia> creator("PANPHASIA");
|
||||
|
|
|
|||
|
|
@ -1,189 +0,0 @@
|
|||
// This file is part of monofonIC (MUSIC2)
|
||||
// A software package to generate ICs for cosmological simulations
|
||||
// Copyright (C) 2021 by Oliver Hahn and Adrian Jenkins (this file)
|
||||
// but see distinct licensing for PANPHASIA below
|
||||
//
|
||||
// monofonIC is free software: you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation, either version 3 of the License, or
|
||||
// (at your option) any later version.
|
||||
//
|
||||
// monofonIC is distributed in the hope that it will be useful,
|
||||
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
// GNU General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
//
|
||||
// IMPORTANT NOTICE:
|
||||
// Note that PANPHASIA itself is not released under the GPL. Make sure
|
||||
// to read and agree to its distinct licensing before you use or modify
|
||||
// the code below or in the /external/panphasia directory which can be
|
||||
// found here: http://icc.dur.ac.uk/Panphasia.php
|
||||
// NOTE THAT PANPHASIA REQUIRES REGISTRATION ON THIS WEBSITE PRIOR TO USE
|
||||
|
||||
#if defined(USE_PANPHASIA_HO)
|
||||
|
||||
#include <general.hh>
|
||||
#include <random_plugin.hh>
|
||||
#include <config_file.hh>
|
||||
|
||||
#include <vector>
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
|
||||
#ifdef _OPENMP
|
||||
#include <omp.h>
|
||||
#endif
|
||||
|
||||
#include <grid_fft.hh>
|
||||
|
||||
namespace PANPHASIA2
|
||||
{
|
||||
extern "C"
|
||||
{
|
||||
#include <panphasia_ho/panphasia_functions.h>
|
||||
extern size_t descriptor_base_size;
|
||||
}
|
||||
}
|
||||
|
||||
#include "PANPHASIA.hh"
|
||||
|
||||
class RNG_panphasia_ho : public RNG_plugin
|
||||
{
|
||||
private:
|
||||
protected:
|
||||
std::string descriptor_string_;
|
||||
int num_threads_;
|
||||
int panphasia_mode_;
|
||||
size_t grid_res_;
|
||||
real_t boxlength_;
|
||||
|
||||
PANPHASIA1::RNG *ppan1_rng_;
|
||||
|
||||
public:
|
||||
explicit RNG_panphasia_ho(config_file &cf) : RNG_plugin(cf)
|
||||
{
|
||||
|
||||
#ifdef _OPENMP
|
||||
num_threads_ = omp_get_max_threads();
|
||||
#else
|
||||
num_threads_ = 1;
|
||||
#endif
|
||||
|
||||
descriptor_string_ = pcf_->get_value<std::string>("random", "descriptor");
|
||||
grid_res_ = pcf_->get_value<size_t>("setup", "GridRes");
|
||||
boxlength_ = pcf_->get_value<real_t>("setup", "BoxLength");
|
||||
|
||||
if( pcf_->get_value<bool>("setup", "DoFixing") ){
|
||||
music::flog << "Fixing all the modes to the mean power negates any advantage of using the Panphasia field.\n";
|
||||
music::flog << "With the panphasia_ho it is possible by choosing the descriptor to fix the largest modes without losing the ability to resimulate to much higher resolution.\n";
|
||||
throw std::runtime_error("PANPHASIA_HO: incompatible parameter.");
|
||||
}
|
||||
|
||||
panphasia_mode_ = 0;
|
||||
PANPHASIA2::parse_and_validate_descriptor_(descriptor_string_.c_str(), &panphasia_mode_);
|
||||
|
||||
if (panphasia_mode_ == 0)
|
||||
{
|
||||
ppan1_rng_ = new PANPHASIA1::RNG(&cf);
|
||||
}
|
||||
}
|
||||
|
||||
~RNG_panphasia_ho()
|
||||
{
|
||||
if (panphasia_mode_ == 0)
|
||||
{
|
||||
delete ppan1_rng_;
|
||||
}
|
||||
}
|
||||
|
||||
bool isMultiscale() const { return true; }
|
||||
|
||||
void Run_Panphasia_Highorder(Grid_FFT<real_t> &g);
|
||||
|
||||
void Fill_Grid(Grid_FFT<real_t> &g)
|
||||
{
|
||||
switch (panphasia_mode_)
|
||||
{
|
||||
|
||||
case 0: // old mode
|
||||
music::ilog << "PANPHASIA: Old descriptor" << std::endl;
|
||||
ppan1_rng_->Fill(g);
|
||||
break;
|
||||
|
||||
case 1: // PANPHASIA HO descriptor
|
||||
music::ilog << "PANPHASIA: New descriptor" << std::endl;
|
||||
this->Run_Panphasia_Highorder(g);
|
||||
break;
|
||||
|
||||
default: // unknown PANPHASIA mode
|
||||
music::elog << "PANPHASIA: Something went wrong with descriptor" << std::endl;
|
||||
abort();
|
||||
break;
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
void RNG_panphasia_ho::Run_Panphasia_Highorder(Grid_FFT<real_t> &g)
|
||||
{
|
||||
int verbose = 0;
|
||||
int error;
|
||||
size_t x0 = 0, y0 = 0, z0 = 0;
|
||||
size_t rel_level;
|
||||
int fdim = 2; //Option to scale Fourier grid dimension relative to Panphasia coefficient grid
|
||||
|
||||
//char descriptor[300] = "[Panph6,L20,(424060,82570,148256),S1,KK0,CH-999,Auriga_100_vol2]";
|
||||
|
||||
PANPHASIA2::PANPHASIA_init_descriptor_(descriptor_string_.c_str(), &verbose);
|
||||
|
||||
printf("Descriptor %s\n ngrid_load %lu\n", descriptor_string_.c_str(), grid_res_);
|
||||
|
||||
// Choose smallest value of level to equal of exceed grid_res_)
|
||||
|
||||
for (rel_level = 0; fdim * (PANPHASIA2::descriptor_base_size <<rel_level) < grid_res_; rel_level++);
|
||||
|
||||
printf("Setting relative level = %lu\n", rel_level);
|
||||
|
||||
if ((error = PANPHASIA2::PANPHASIA_init_level_(&rel_level, &x0, &y0, &z0, &verbose)))
|
||||
{
|
||||
printf("Abort: PANPHASIA_init_level_ :error code %d\n", error);
|
||||
abort();
|
||||
};
|
||||
|
||||
//======================= FFTW ==============================
|
||||
|
||||
ptrdiff_t alloc_local, local_n0, local_0_start;
|
||||
|
||||
size_t N0 = fdim * (PANPHASIA2::descriptor_base_size << rel_level);
|
||||
|
||||
alloc_local = FFTW_MPI_LOCAL_SIZE_3D(N0, N0, N0 + 2, MPI_COMM_WORLD, &local_n0, &local_0_start);
|
||||
|
||||
Grid_FFT<real_t> pan_grid({{N0, N0, N0}}, {{boxlength_, boxlength_, boxlength_}});
|
||||
|
||||
assert(pan_grid.n_[0] == N0);
|
||||
assert(pan_grid.n_[1] == N0);
|
||||
assert(pan_grid.n_[2] == N0);
|
||||
assert(pan_grid.local_0_start_ == local_0_start);
|
||||
assert(pan_grid.local_0_size_ == local_n0);
|
||||
assert(pan_grid.ntot_ == size_t(alloc_local));
|
||||
|
||||
pan_grid.FourierTransformForward(false);
|
||||
|
||||
FFTW_COMPLEX *Panphasia_White_Noise_Field = reinterpret_cast<FFTW_COMPLEX *>(&pan_grid.data_[0]);
|
||||
|
||||
// Panphasia_White_Noise_Field = FFTW_ALLOC_COMPLEX(alloc_local);
|
||||
|
||||
if ((error = PANPHASIA2::PANPHASIA_compute_kspace_field_(rel_level, N0, local_n0, local_0_start, Panphasia_White_Noise_Field)))
|
||||
{
|
||||
music::elog << "Error code from PANPHASIA_compute ... (ErrCode = " << error << ")" << std::endl;
|
||||
};
|
||||
|
||||
pan_grid.FourierInterpolateCopyTo(g);
|
||||
}
|
||||
namespace
|
||||
{
|
||||
RNG_plugin_creator_concrete<RNG_panphasia_ho> creator("PANPHASIA_HO");
|
||||
}
|
||||
#endif // defined(USE_PANPHASIA_HO)
|
||||
Loading…
Reference in a new issue