added option to normalise mismatch in mean density upward or downard or not. defaults to not for panphasia now

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<real_t> &u, bool kspace)
 {
-	unsigned levelmin_TF = u.levelmin();
+	const unsigned levelmin_TF = u.levelmin();
-	bool benforce_coarse = true;//!kspace;
+	const bool benforce_coarse = !kspace;
 
 	if (kspace)
 	{
@@ -509,4 +549,7 @@ void coarsen_density(const refinement_hierarchy &rh, GridHierarchy<real_t> &u, b
 						rh.size(i, 0), rh.size(i, 1), rh.size(i, 2), benforce_coarse);
 		}
 	}
+	if( !benforce_coarse ){
+		normalize_levelmin_density( u );
+	}
 }

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"))
+	const bool use_fourier_coarsening = cf.get_value_safe<bool>("setup", "fourier_splicing", true);
-		cf.insert_value("setup", "kspace_TF", "yes");
 
-	bool bspectral_sampling = true;//cf.get_value<bool>("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);
 

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);
 			}
 		}

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() {}