From c184accf00db04d3a729644bfd93fe8e72f41664 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Mon, 13 Feb 2023 15:07:54 -0800 Subject: [PATCH] WIP panphasia field now agrees within some errors with the previous result. --- src/convolution_kernel.cc | 3 +- src/densities.cc | 14 ++++-- src/mesh.hh | 89 +++++++++++++++++++++++++-------- src/plugins/random_music.cc | 2 + src/plugins/random_panphasia.cc | 29 ++++++++--- 5 files changed, 103 insertions(+), 34 deletions(-) diff --git a/src/convolution_kernel.cc b/src/convolution_kernel.cc index 3aa1f3d..c766e37 100644 --- a/src/convolution_kernel.cc +++ b/src/convolution_kernel.cc @@ -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(pd); diff --git a/src/densities.cc b/src/densities.cc index 589fe33..5ce7a35 100644 --- a/src/densities.cc +++ b/src/densities.cc @@ -445,6 +445,7 @@ void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type t bool fix = cf.getValueSafe("setup","fix_mode_amplitude",false); bool flip = cf.getValueSafe("setup","flip_mode_amplitude",false); + bool fourier_splicing = cf.getValueSafe("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(the_tf_kernel->fetch_kernel(levelmin + i, true), reinterpret_cast(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); } diff --git a/src/mesh.hh b/src/mesh.hh index efb82a5..744be4d 100644 --- a/src/mesh.hh +++ b/src/mesh.hh @@ -1208,11 +1208,12 @@ class refinement_hierarchy std::vector 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<( 1, (1<<(levelmin+1)) / (2*gcd(U,1<<(levelmin+1) )) );*/ + + int level_m = 0; + while( base_unit * (1<( 1, (1<("setup", "no_shift", false); bool force_shift = cf_.getValueSafe("setup", "force_shift", false); + gridding_unit_ = cf.getValueSafe("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("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) diff --git a/src/plugins/random_music.cc b/src/plugins/random_music.cc index 0bfdf6b..035288e 100644 --- a/src/plugins/random_music.cc +++ b/src/plugins/random_music.cc @@ -50,6 +50,8 @@ public: disk_cached_ = pcf_->getValueSafe("random", "disk_cached", true); restart_ = pcf_->getValueSafe("random", "restart", false); + pcf_->insertValue("setup","fourier_splicing","true"); + mem_cache_.assign(levelmax_ - levelmin_ + 1, (std::vector *)NULL); if (restart_ && !disk_cached_) diff --git a/src/plugins/random_panphasia.cc b/src/plugins/random_panphasia.cc index 53dc456..0fa08fd 100644 --- a/src/plugins/random_panphasia.cc +++ b/src/plugins/random_panphasia.cc @@ -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 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("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 &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 &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 &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 &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 &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 &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<