mirror of
https://github.com/cosmo-sims/MUSIC.git
synced 2024-09-19 17:03:46 +02:00
WIP panphasia field now agrees within some errors with the previous result.
This commit is contained in:
parent
ced852069e
commit
c184accf00
5 changed files with 103 additions and 34 deletions
|
@ -34,7 +34,8 @@ void perform(kernel *pk, void *pd, bool shift, bool fix, bool flip)
|
|||
double fftnormp = 1.0/sqrt((double)cparam_.nx * (double)cparam_.ny * (double)cparam_.nz);
|
||||
double fftnorm = pow(2.0 * M_PI, 1.5) / sqrt(cparam_.lx * cparam_.ly * cparam_.lz) * fftnormp;
|
||||
|
||||
fftw_complex *cdata, *ckernel;
|
||||
fftw_complex *cdata;
|
||||
[[maybe_unused]] fftw_complex *ckernel;
|
||||
fftw_real *data;
|
||||
|
||||
data = reinterpret_cast<fftw_real *>(pd);
|
||||
|
|
|
@ -445,6 +445,7 @@ void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type t
|
|||
|
||||
bool fix = cf.getValueSafe<bool>("setup","fix_mode_amplitude",false);
|
||||
bool flip = cf.getValueSafe<bool>("setup","flip_mode_amplitude",false);
|
||||
bool fourier_splicing = cf.getValueSafe<bool>("setup","fourier_splicing",true);
|
||||
|
||||
if( fix && levelmin != levelmax ){
|
||||
LOGWARN("You have chosen mode fixing for a zoom. This is not well tested,\n please proceed at your own risk...");
|
||||
|
@ -511,10 +512,12 @@ void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type t
|
|||
convolution::perform<real_t>(the_tf_kernel->fetch_kernel(levelmin + i, true),
|
||||
reinterpret_cast<void *>(fine->get_data_ptr()), shift, fix, flip);
|
||||
|
||||
if (i == 1)
|
||||
fft_interpolate(*top, *fine, true);
|
||||
else
|
||||
fft_interpolate(*coarse, *fine, false);
|
||||
if( fourier_splicing ){
|
||||
if (i == 1)
|
||||
fft_interpolate(*top, *fine, true);
|
||||
else
|
||||
fft_interpolate(*coarse, *fine, false);
|
||||
}
|
||||
|
||||
delta.add_patch(refh.offset(levelmin + i, 0),
|
||||
refh.offset(levelmin + i, 1),
|
||||
|
@ -548,6 +551,9 @@ void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type t
|
|||
std::cout << " - Density calculation took " << tend - tstart << "s." << std::endl;
|
||||
#endif
|
||||
|
||||
if( !fourier_splicing ){
|
||||
coarsen_density(refh,delta,false);
|
||||
}
|
||||
LOGUSER("Finished computing the density field in %fs", tend - tstart);
|
||||
}
|
||||
|
||||
|
|
89
src/mesh.hh
89
src/mesh.hh
|
@ -1208,11 +1208,12 @@ class refinement_hierarchy
|
|||
std::vector<index3_t> len_;
|
||||
|
||||
unsigned
|
||||
levelmin_, //!< minimum grid level for Poisson solver
|
||||
levelmax_, //!< maximum grid level for all operations
|
||||
levelmin_tf_, //!< minimum grid level for density calculation
|
||||
padding_, //!< padding in number of coarse cells between refinement levels
|
||||
blocking_factor_;
|
||||
levelmin_, //!< minimum grid level for Poisson solver
|
||||
levelmax_, //!< maximum grid level for all operations
|
||||
levelmin_tf_, //!< minimum grid level for density calculation
|
||||
padding_, //!< padding in number of coarse cells between refinement levels
|
||||
blocking_factor_, //!< blocking factor of grids, necessary fo BoxLib codes such as NyX
|
||||
gridding_unit_; //!< internal blocking factor of grids, necessary for Panphasia
|
||||
|
||||
int margin_; //!< number of cells used for additional padding for convolutions with isolated boundaries (-1 = double padding)
|
||||
|
||||
|
@ -1231,6 +1232,30 @@ class refinement_hierarchy
|
|||
index3_t xshift_; //!< shift of refinement region in coarse cells (in order to center it in the domain)
|
||||
double rshift_[3];
|
||||
|
||||
//! calculates the greatest common divisor
|
||||
int gcd(int a, int b) const {
|
||||
return b == 0 ? a : gcd(b, a % b);
|
||||
}
|
||||
|
||||
// calculates the cell shift in units of levelmin grid cells if there is an additional constraint to be
|
||||
// congruent with another grid that partitions the same space in multiples of "base_unit"
|
||||
int get_shift_unit( int base_unit, int levelmin ) const {
|
||||
/*int Lp = 0;
|
||||
while( base_unit * (1<<Lp) < 1<<(levelmin+1) ){
|
||||
++Lp;
|
||||
}
|
||||
int U = base_unit * (1<<Lp);
|
||||
|
||||
return std::max<int>( 1, (1<<(levelmin+1)) / (2*gcd(U,1<<(levelmin+1) )) );*/
|
||||
|
||||
int level_m = 0;
|
||||
while( base_unit * (1<<level_m) < (1<<levelmin) )
|
||||
++level_m;
|
||||
|
||||
return std::max<int>( 1, (1<<levelmin)/gcd(base_unit * (1<<level_m),(1<<levelmin)) );
|
||||
|
||||
}
|
||||
|
||||
public:
|
||||
//! copy constructor
|
||||
refinement_hierarchy(const refinement_hierarchy &rh)
|
||||
|
@ -1256,6 +1281,16 @@ public:
|
|||
bool bnoshift = cf_.getValueSafe<bool>("setup", "no_shift", false);
|
||||
bool force_shift = cf_.getValueSafe<bool>("setup", "force_shift", false);
|
||||
|
||||
gridding_unit_ = cf.getValueSafe<unsigned>("setup", "gridding_unit", 2);
|
||||
|
||||
if (gridding_unit_ != 2 && blocking_factor_==0) {
|
||||
blocking_factor_ = gridding_unit_; // THIS WILL LIKELY CAUSE PROBLEMS WITH NYX
|
||||
}else if (gridding_unit_ != 2 && blocking_factor_!=0 && gridding_unit_!=blocking_factor_ ) {
|
||||
LOGERR("incompatible gridding unit %d and blocking factor specified", gridding_unit_, blocking_factor_ );
|
||||
throw std::runtime_error("Incompatible gridding unit and blocking factor!");
|
||||
}
|
||||
|
||||
|
||||
//... call the region generator
|
||||
if (levelmin_ != levelmax_)
|
||||
{
|
||||
|
@ -1292,9 +1327,18 @@ public:
|
|||
|
||||
if ((levelmin_ != levelmax_) && (!bnoshift || force_shift))
|
||||
{
|
||||
xshift_[0] = (int)((0.5 - xc[0]) * ncoarse);
|
||||
xshift_[1] = (int)((0.5 - xc[1]) * ncoarse);
|
||||
xshift_[2] = (int)((0.5 - xc[2]) * ncoarse);
|
||||
int random_base_grid_unit = cf.getValueSafe<int>("random","base_unit",1);
|
||||
int shift_unit = get_shift_unit( random_base_grid_unit, levelmin_ );
|
||||
if( shift_unit != 1 ){
|
||||
LOGINFO("volume can only be shifted by multiples of %d coarse cells.",shift_unit);
|
||||
}
|
||||
xshift_[0] = (int)((0.5-xc[0]) * (double)ncoarse / shift_unit + 0.5) * shift_unit;//ARJ(int)((0.5 - xc[0]) * ncoarse);
|
||||
xshift_[1] = (int)((0.5-xc[1]) * (double)ncoarse / shift_unit + 0.5) * shift_unit;//ARJ(int)((0.5 - xc[1]) * ncoarse);
|
||||
xshift_[2] = (int)((0.5-xc[2]) * (double)ncoarse / shift_unit + 0.5) * shift_unit;//ARJ(int)((0.5 - xc[2]) * ncoarse);
|
||||
|
||||
// xshift_[0] = (int)((0.5 - xc[0]) * ncoarse);
|
||||
// xshift_[1] = (int)((0.5 - xc[1]) * ncoarse);
|
||||
// xshift_[2] = (int)((0.5 - xc[2]) * ncoarse);
|
||||
}
|
||||
else
|
||||
{
|
||||
|
@ -1414,12 +1458,15 @@ public:
|
|||
else
|
||||
{
|
||||
//... require alignment with coarser grid
|
||||
il -= il % 2;
|
||||
jl -= jl % 2;
|
||||
kl -= kl % 2;
|
||||
ir += ir % 2;
|
||||
jr += jr % 2;
|
||||
kr += kr % 2;
|
||||
LOGINFO("Internal refinement bounding box error: [%d,%d]x[%d,%d]x[%d,%d]", il, ir, jl, jr, kl, kr);
|
||||
|
||||
il -= il % gridding_unit_;
|
||||
jl -= jl % gridding_unit_;
|
||||
kl -= kl % gridding_unit_;
|
||||
|
||||
ir = ((ir%gridding_unit_)!=0)? (ir/gridding_unit_ + 1)*gridding_unit_ : ir;
|
||||
jr = ((jr%gridding_unit_)!=0)? (jr/gridding_unit_ + 1)*gridding_unit_ : jr;
|
||||
kr = ((kr%gridding_unit_)!=0)? (kr/gridding_unit_ + 1)*gridding_unit_ : kr;
|
||||
}
|
||||
|
||||
// if doing unigrid, set region to whole box
|
||||
|
@ -1511,7 +1558,6 @@ public:
|
|||
jr = (int)((double)jr / nref + 1.0) * nref;
|
||||
kr = (int)((double)kr / nref + 1.0) * nref;
|
||||
}
|
||||
|
||||
else if (preserve_dims_)
|
||||
{
|
||||
//... require alignment with coarser grid
|
||||
|
@ -1528,12 +1574,13 @@ public:
|
|||
else
|
||||
{
|
||||
//... require alignment with coarser grid
|
||||
il -= il % 2;
|
||||
jl -= jl % 2;
|
||||
kl -= kl % 2;
|
||||
ir += ir % 2;
|
||||
jr += jr % 2;
|
||||
kr += kr % 2;
|
||||
il -= il % gridding_unit_;
|
||||
jl -= jl % gridding_unit_;
|
||||
kl -= kl % gridding_unit_;
|
||||
|
||||
ir = ((ir%gridding_unit_)!=0)? (ir/gridding_unit_ + 1)*gridding_unit_ : ir;
|
||||
jr = ((jr%gridding_unit_)!=0)? (jr/gridding_unit_ + 1)*gridding_unit_ : jr;
|
||||
kr = ((kr%gridding_unit_)!=0)? (kr/gridding_unit_ + 1)*gridding_unit_ : kr;
|
||||
}
|
||||
|
||||
if (il >= ir || jl >= jr || kl >= kr || il < 0 || jl < 0 || kl < 0)
|
||||
|
|
|
@ -50,6 +50,8 @@ public:
|
|||
disk_cached_ = pcf_->getValueSafe<bool>("random", "disk_cached", true);
|
||||
restart_ = pcf_->getValueSafe<bool>("random", "restart", false);
|
||||
|
||||
pcf_->insertValue("setup","fourier_splicing","true");
|
||||
|
||||
mem_cache_.assign(levelmax_ - levelmin_ + 1, (std::vector<real_t> *)NULL);
|
||||
|
||||
if (restart_ && !disk_cached_)
|
||||
|
|
|
@ -92,6 +92,7 @@ protected:
|
|||
int coordinate_system_shift_[3];
|
||||
int ix_abs_[3], ix_per_[3], ix_rel_[3], level_p_, lextra_;
|
||||
const refinement_hierarchy *prefh_;
|
||||
std::array<int,3> margins_;
|
||||
|
||||
struct panphasia_descriptor
|
||||
{
|
||||
|
@ -180,6 +181,8 @@ public:
|
|||
ss.str(std::string());
|
||||
ss << pdescriptor_->i_base;
|
||||
pcf_->insertValue("random", "base_unit", ss.str());
|
||||
|
||||
pcf_->insertValue("setup","fourier_splicing","false");
|
||||
}
|
||||
|
||||
void initialize_for_grid_structure(const refinement_hierarchy &refh)
|
||||
|
@ -189,6 +192,12 @@ public:
|
|||
levelmin_final_ = pcf_->getValue<unsigned>("setup", "levelmin");
|
||||
levelmax_ = prefh_->levelmax();
|
||||
|
||||
if( refh.get_margin() < 0 ){
|
||||
margins_ = {-1,-1,-1};
|
||||
}else{
|
||||
margins_ = { refh.get_margin(), refh.get_margin(), refh.get_margin() };
|
||||
}
|
||||
|
||||
clear_panphasia_thread_states();
|
||||
LOGINFO("PANPHASIA: running with %d threads", num_threads_);
|
||||
|
||||
|
@ -318,7 +327,13 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
|
|||
}
|
||||
else
|
||||
{
|
||||
ileft_corner[k] = (ileft[k] - nx[k] / 4 + (1 << level)) % (1 << level); // Isolated
|
||||
if( margins_[0] < 0 ){
|
||||
ileft_corner[k] = (ileft[k] - nx[k] / 4 + (1 << level)) % (1 << level); // Isolated
|
||||
ileft_corner[k] = (ileft[k] - nx[k] / 4 + (1 << level)) % (1 << level); // Isolated
|
||||
}else{
|
||||
ileft_corner[k] = (ileft[k] - margins_[k] + (1 << level)) % (1 << level); // Isolated
|
||||
ileft_corner[k] = (ileft[k] - margins_[k] + (1 << level)) % (1 << level); // Isolated
|
||||
}
|
||||
}
|
||||
iexpand_left[k] = (ileft_corner[k] % grid_m_ == 0) ? 0 : ileft_corner[k] % grid_m_;
|
||||
// fprintf(stderr, "dim=%c : ileft = %d, ileft_corner %d, nx = %d\n", 'x' + k, ileft[k],ileft_corner[k],nx[k]);
|
||||
|
@ -349,7 +364,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
|
|||
LOGERR("Fatal error: non-cubic refinement being requested");
|
||||
|
||||
inter_grid_phase_adjustment_ = M_PI * (1.0 / (double)nx_m[0] - 1.0 / (double)nxremap[0]);
|
||||
LOGUSER("The value of the phase adjustement is %f\n", inter_grid_phase_adjustment_);
|
||||
// LOGINFO("The value of the phase adjustement is %f\n", inter_grid_phase_adjustment_);
|
||||
|
||||
// LOGINFO("ileft[0],ileft[1],ileft[2] %d %d %d", ileft[0], ileft[1], ileft[2]);
|
||||
// LOGINFO("ileft_corner[0,1,2] %d %d %d", ileft_corner[0], ileft_corner[1], ileft_corner[2]);
|
||||
|
@ -425,8 +440,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
|
|||
set_phases_and_rel_origin_(&lstate[mythread], descriptor, &level_p, &ix_rel[0], &ix_rel[1], &ix_rel[2],
|
||||
&verbosity);
|
||||
|
||||
LOGUSER(" called set_phases_and_rel_origin level %d ix_rel iy_rel iz_rel %d %d %d\n", level_p, ix_rel[0],
|
||||
ix_rel[1], ix_rel[2]);
|
||||
// LOGUSER(" called set_phases_and_rel_origin level %d ix_rel iy_rel iz_rel %d %d %d\n", level_p, ix_rel[0], ix_rel[1], ix_rel[2]);
|
||||
|
||||
odd_x = ix_rel[0] % 2;
|
||||
odd_y = ix_rel[1] % 2;
|
||||
|
@ -620,8 +634,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
|
|||
set_phases_and_rel_origin_(&lstate[mythread], descriptor, &level_p, &ix_rel[0], &ix_rel[1], &ix_rel[2],
|
||||
&verbosity);
|
||||
|
||||
LOGUSER(" called set_phases_and_rel_origin level %d ix_rel iy_rel iz_rel %d %d %d\n", level_p, ix_rel[0],
|
||||
ix_rel[1], ix_rel[2]);
|
||||
// LOGINFO(" called set_phases_and_rel_origin level %d ix_rel iy_rel iz_rel %d %d %d\n", level_p, ix_rel[0], ix_rel[1], ix_rel[2]);
|
||||
|
||||
odd_x = ix_rel[0] % 2;
|
||||
odd_y = ix_rel[1] % 2;
|
||||
|
@ -769,7 +782,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
|
|||
if (incongruent_fields_)
|
||||
{
|
||||
|
||||
LOGINFO("Remapping fields from dimension %d -> %d", nxremap[0], nx_m[0]);
|
||||
LOGUSER("Remapping fields from dimension %d -> %d", nxremap[0], nx_m[0]);
|
||||
memset(pr1, 0, ngp * sizeof(fftw_real));
|
||||
|
||||
#pragma omp parallel for
|
||||
|
@ -824,7 +837,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
|
|||
delete[] pr3;
|
||||
delete[] pr4;
|
||||
|
||||
LOGINFO("Copying random field data %d,%d,%d -> %d,%d,%d", nxremap[0], nxremap[1], nxremap[2], nx[0], nx[1], nx[2]);
|
||||
LOGUSER("Copying random field data %d,%d,%d -> %d,%d,%d", nxremap[0], nxremap[1], nxremap[2], nx[0], nx[1], nx[2]);
|
||||
|
||||
// n = 1<<level;
|
||||
// ng = n;
|
||||
|
|
Loading…
Reference in a new issue