From 74a1727d3d4f6e77aef79810904c16e4276e9905 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Wed, 8 Mar 2023 20:07:36 +0100 Subject: [PATCH] added option to normalise mismatch in mean density upward or downard or not. defaults to not for panphasia now --- src/densities.cc | 47 +++++++++++++++++++++++++++++++++++-- src/main.cc | 29 +++++++++-------------- src/mesh.hh | 28 ++++++++++++++++++++++ src/plugins/random_music.cc | 5 +++- 4 files changed, 88 insertions(+), 21 deletions(-) diff --git a/src/densities.cc b/src/densities.cc index 10b4ecc..ad51991 100644 --- a/src/densities.cc +++ b/src/densities.cc @@ -484,10 +484,50 @@ void normalize_density(grid_hierarchy &delta) } } +void normalize_levelmin_density(grid_hierarchy &delta) +{ + //return; + + long double sum = 0.0; + unsigned levelmin = delta.levelmin(), levelmax = delta.levelmax(); + + { + size_t nx, ny, nz; + + nx = delta.get_grid(levelmin)->size(0); + ny = delta.get_grid(levelmin)->size(1); + nz = delta.get_grid(levelmin)->size(2); + +#pragma omp parallel for reduction(+ : sum) + for (int ix = 0; ix < (int)nx; ++ix) + for (size_t iy = 0; iy < ny; ++iy) + for (size_t iz = 0; iz < nz; ++iz) + sum += (*delta.get_grid(levelmin))(ix, iy, iz); + + sum /= (double)(nx * ny * nz); + } + + music::ilog << "- Top grid mean density is off by " << sum << ", correcting..." << std::endl; + + { + size_t nx, ny, nz; + nx = delta.get_grid(levelmin)->size(0); + ny = delta.get_grid(levelmin)->size(1); + nz = delta.get_grid(levelmin)->size(2); + +#pragma omp parallel for + for (int ix = 0; ix < (int)nx; ++ix) + for (size_t iy = 0; iy < ny; ++iy) + for (size_t iz = 0; iz < nz; ++iz) + (*delta.get_grid(levelmin))(ix, iy, iz) -= sum; + } +} + + void coarsen_density(const refinement_hierarchy &rh, GridHierarchy &u, bool kspace) { - unsigned levelmin_TF = u.levelmin(); - bool benforce_coarse = true;//!kspace; + const unsigned levelmin_TF = u.levelmin(); + const bool benforce_coarse = !kspace; if (kspace) { @@ -509,4 +549,7 @@ void coarsen_density(const refinement_hierarchy &rh, GridHierarchy &u, b rh.size(i, 0), rh.size(i, 1), rh.size(i, 2), benforce_coarse); } } + if( !benforce_coarse ){ + normalize_levelmin_density( u ); + } } diff --git a/src/main.cc b/src/main.cc index 629b000..0c8a32f 100644 --- a/src/main.cc +++ b/src/main.cc @@ -453,15 +453,8 @@ int main(int argc, const char *argv[]) } // .. determine if spectral sampling should be used - if (!cf.contains_key("setup", "kspace_TF")) - cf.insert_value("setup", "kspace_TF", "yes"); + const bool use_fourier_coarsening = cf.get_value_safe("setup", "fourier_splicing", true); - bool bspectral_sampling = true;//cf.get_value("setup", "kspace_TF"); - - if (bspectral_sampling) - music::ilog.Print("Using k-space sampled transfer functions..."); - else - music::ilog.Print("Using real space sampled transfer functions..."); //------------------------------------------------------------------------------ //... initialize cosmology @@ -622,7 +615,7 @@ int main(int argc, const char *argv[]) my_tf_type = delta_matter; GenerateDensityHierarchy(cf, the_cosmo_calc.get(), my_tf_type, rh_TF, rand, f, false, false); - coarsen_density(rh_Poisson, f, bspectral_sampling); + coarsen_density(rh_Poisson, f, use_fourier_coarsening); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -686,7 +679,7 @@ int main(int argc, const char *argv[]) music::ilog << "-------------------------------------------------------------------------------" << std::endl; music::ulog.Print("Computing baryon density..."); GenerateDensityHierarchy(cf, the_cosmo_calc.get(), delta_baryon, rh_TF, rand, f, false, bbshift); - coarsen_density(rh_Poisson, f, bspectral_sampling); + coarsen_density(rh_Poisson, f, use_fourier_coarsening); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -759,7 +752,7 @@ int main(int argc, const char *argv[]) { music::ulog.Print("Generating velocity perturbations..."); GenerateDensityHierarchy(cf, the_cosmo_calc.get(), theta_cdm, rh_TF, rand, f, false, false); - coarsen_density(rh_Poisson, f, bspectral_sampling); + coarsen_density(rh_Poisson, f, use_fourier_coarsening); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); u = f; @@ -829,7 +822,7 @@ int main(int argc, const char *argv[]) //... we do baryons and have velocity transfer functions, or we do SPH and not to shift //... do DM first GenerateDensityHierarchy(cf, the_cosmo_calc.get(), theta_cdm, rh_TF, rand, f, false, false); - coarsen_density(rh_Poisson, f, bspectral_sampling); + coarsen_density(rh_Poisson, f, use_fourier_coarsening); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -888,7 +881,7 @@ int main(int argc, const char *argv[]) music::ulog.Print("Computing baryon velocitites..."); //... do baryons GenerateDensityHierarchy(cf, the_cosmo_calc.get(), theta_baryon, rh_TF, rand, f, false, bbshift); - coarsen_density(rh_Poisson, f, bspectral_sampling); + coarsen_density(rh_Poisson, f, use_fourier_coarsening); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -970,7 +963,7 @@ int main(int argc, const char *argv[]) music::ilog << "-------------------------------------------------------------------------------" << std::endl; GenerateDensityHierarchy(cf, the_cosmo_calc.get(), my_tf_type, rh_TF, rand, f, false, false); - coarsen_density(rh_Poisson, f, bspectral_sampling); + coarsen_density(rh_Poisson, f, use_fourier_coarsening); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -1076,7 +1069,7 @@ int main(int argc, const char *argv[]) music::ulog.Print("Computing baryon displacements..."); GenerateDensityHierarchy(cf, the_cosmo_calc.get(), theta_baryon, rh_TF, rand, f, false, bbshift); - coarsen_density(rh_Poisson, f, bspectral_sampling); + coarsen_density(rh_Poisson, f, use_fourier_coarsening); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -1174,7 +1167,7 @@ int main(int argc, const char *argv[]) my_tf_type = delta_matter; GenerateDensityHierarchy(cf, the_cosmo_calc.get(), my_tf_type, rh_TF, rand, f, false, false); - coarsen_density(rh_Poisson, f, bspectral_sampling); + coarsen_density(rh_Poisson, f, use_fourier_coarsening); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -1272,7 +1265,7 @@ int main(int argc, const char *argv[]) music::ulog.Print("Computing baryon density..."); GenerateDensityHierarchy(cf, the_cosmo_calc.get(), delta_baryon, rh_TF, rand, f, true, false); - coarsen_density(rh_Poisson, f, bspectral_sampling); + coarsen_density(rh_Poisson, f, use_fourier_coarsening); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -1315,7 +1308,7 @@ int main(int argc, const char *argv[]) music::ulog.Print("Computing baryon displacements..."); GenerateDensityHierarchy(cf, the_cosmo_calc.get(), delta_baryon, rh_TF, rand, f, false, bbshift); - coarsen_density(rh_Poisson, f, bspectral_sampling); + coarsen_density(rh_Poisson, f, use_fourier_coarsening); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); diff --git a/src/mesh.hh b/src/mesh.hh index af43f90..11f52f6 100644 --- a/src/mesh.hh +++ b/src/mesh.hh @@ -1169,6 +1169,34 @@ public: for (unsigned k = 0; k < nz; ++k) (*mnew)(i, j, k) += (coarsesum - finesum); + music::ilog.Print(" .level %d : corrected patch overlap mean value by %f", ilevel, coarsesum - finesum); + } + }else{ + //... enforce fine mean density over same patch + if (ilevel > levelmin()) + { + int ox = m_pgrids[ilevel]->offset(0); + int oy = m_pgrids[ilevel]->offset(1); + int oz = m_pgrids[ilevel]->offset(2); + + #pragma omp parallel for reduction(+:coarsesum,coarsecount) collapse(3) + for (unsigned i = 0; i < nx / 2; ++i) + for (unsigned j = 0; j < ny / 2; ++j) + for (unsigned k = 0; k < nz / 2; ++k) + { + coarsesum += (*m_pgrids[ilevel - 1])(i + ox, j + oy, k + oz); + coarsecount++; + } + + coarsesum /= (double)coarsecount; + finesum /= (double)finecount; + + #pragma omp parallel for collapse(3) + for (unsigned i = 0; i < nx / 2; ++i) + for (unsigned j = 0; j < ny / 2; ++j) + for (unsigned k = 0; k < nz / 2; ++k) + (*m_pgrids[ilevel - 1])(i + ox, j + oy, k + oz) -= (coarsesum - finesum); + music::ilog.Print(" .level %d : corrected patch overlap mean value by %f", ilevel, coarsesum - finesum); } } diff --git a/src/plugins/random_music.cc b/src/plugins/random_music.cc index d935258..af35b3e 100644 --- a/src/plugins/random_music.cc +++ b/src/plugins/random_music.cc @@ -34,7 +34,10 @@ protected: void store_rnd(int ilevel, rng *prng); public: - explicit RNG_music(config_file &cf) : RNG_plugin(cf), initialized_(false) {} + explicit RNG_music(config_file &cf) : RNG_plugin(cf), initialized_(false) + { + pcf_->insert_value("setup","fourier_splicing","true"); + } ~RNG_music() {}