WIP some fixes to noise generators

src/densities.cc

@@ -23,7 +23,7 @@ double Meyer_scaling_function( double k, double kmax )
 	constexpr double fourpithirds{4.0*M_PI/3.0};
 	auto nu = []( double x ){ return x<0.0?0.0:(x<1.0?x:1.0); };
 
-	k = k/kmax * fourpithirds;
+	k = std::abs(k)/kmax * fourpithirds;
 
 	if( k < twopithirds ) return 1.0;
 	else if( k< fourpithirds ){
@@ -126,15 +126,15 @@ void fft_interpolate(m1 &V, m2 &v, bool from_basegrid = false)
 
 	if (!from_basegrid)
 	{
-		oxf +=  mxf - mxf/2;
+		oxf +=  mxf/2;
-		oyf +=  myf - myf/2;
+		oyf +=  myf/2;
-		ozf +=  mzf - mzf/2;
+		ozf +=  mzf/2;
 	}
 	else
 	{
-		oxf -= mxf/2; //nxf / 8; 
+		oxf -= mxf/2;
-		oyf -= myf/2; //nyf / 8; 
+		oyf -= myf/2;
-		ozf -= mzf/2; //nzf / 8; 
+		ozf -= mzf/2;
 	}
 
 	music::ulog.Print("FFT interpolate: offset=%d,%d,%d size=%d,%d,%d", oxf, oyf, ozf, nxf, nyf, nzf);
@@ -185,10 +185,9 @@ void fft_interpolate(m1 &V, m2 &v, bool from_basegrid = false)
 	/*************************************************/
 	//.. perform actual interpolation
 	double fftnorm = 1.0 / ((double)nxf * (double)nyf * (double)nzf);
-	double sqrt8 = 8.0; //sqrt(8.0);
-	double phasefac = -0.5;
 
 	// this enables filtered splicing of coarse and fine modes
+	#pragma omp parallel for
 	for (int i = 0; i < (int)nxc; i++)
 		for (int j = 0; j < (int)nyc; j++)
 			for (int k = 0; k < (int)nzc / 2 + 1; k++)
@@ -210,20 +209,28 @@ void fft_interpolate(m1 &V, m2 &v, bool from_basegrid = false)
 				double ky = (j <= (int)nyc / 2) ? (double)j : (double)(j - (int)nyc);
 				double kz = (k <= (int)nzc / 2) ? (double)k : (double)(k - (int)nzc);
 
-				double phase = phasefac * (kx / nxc + ky / nyc + kz / nzc) * M_PI;
+				double phase = -0.5 * M_PI * (kx / nxc + ky / nyc + kz / nzc);
 
 				std::complex<double> val_phas(cos(phase), sin(phase));
 
 				std::complex<double> val(RE(ccoarse[qc]), IM(ccoarse[qc]));
-				val *= sqrt8 * val_phas;
+				val *= val_phas * 8.0;
+
+				if (false){ //(i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2){
+					double blend_coarse_x = Meyer_scaling_function(kx, nxc / 2);
+					double blend_coarse_y = Meyer_scaling_function(ky, nyc / 2);
+					double blend_coarse_z = Meyer_scaling_function(kz, nzc / 2);
+					double blend_coarse = blend_coarse_x*blend_coarse_y*blend_coarse_z;
+					double blend_fine = 1.0-blend_coarse;
 
-				if (i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2){
-					double blend_coarse = Meyer_scaling_function(kx, nxc / 2) * Meyer_scaling_function(ky, nyc / 2) * Meyer_scaling_function(kz, nzc / 2);
-					double blend_fine = 1.0 - blend_coarse; 
-				
 					RE(cfine[qf]) = blend_fine * RE(cfine[qf]) + blend_coarse * val.real();
 					IM(cfine[qf]) = blend_fine * IM(cfine[qf]) + blend_coarse * val.imag();
 				}
+
+				if (i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2){
+					RE(cfine[qf]) = val.real();
+					IM(cfine[qf]) = val.imag();
+				}
 			}
 
 	delete[] rcoarse;
@@ -451,8 +458,7 @@ void normalize_density(grid_hierarchy &delta)
 		ny = delta.get_grid(levelmin)->size(1);
 		nz = delta.get_grid(levelmin)->size(2);
 
-#pragma omp parallel for reduction(+ \
+#pragma omp parallel for reduction(+ : sum)
-								   : 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)
@@ -481,9 +487,7 @@ void normalize_density(grid_hierarchy &delta)
 void coarsen_density(const refinement_hierarchy &rh, GridHierarchy<real_t> &u, bool kspace)
 {
 	unsigned levelmin_TF = u.levelmin();
-
+	bool benforce_coarse = true;//!kspace;
-	/*for( int i=rh.levelmax(); i>0; --i )
-        mg_straight().restrict( *(u.get_grid(i)), *(u.get_grid(i-1)) );*/
 
 	if (kspace)
 	{
@@ -498,13 +502,11 @@ void coarsen_density(const refinement_hierarchy &rh, GridHierarchy<real_t> &u, b
 
 	for (unsigned i = 1; i <= rh.levelmax(); ++i)
 	{
-		if (rh.offset(i, 0) != u.get_grid(i)->offset(0) || rh.offset(i, 1) != u.get_grid(i)->offset(1) || rh.offset(i, 2) != u.get_grid(i)->offset(2) || rh.size(i, 0) != u.get_grid(i)->size(0) || rh.size(i, 1) != u.get_grid(i)->size(1) || rh.size(i, 2) != u.get_grid(i)->size(2))
+		if (rh.offset(i, 0) != u.get_grid(i)->offset(0) || rh.offset(i, 1) != u.get_grid(i)->offset(1) || rh.offset(i, 2) != u.get_grid(i)->offset(2) 
+			|| rh.size(i, 0) != u.get_grid(i)->size(0) || rh.size(i, 1) != u.get_grid(i)->size(1) || rh.size(i, 2) != u.get_grid(i)->size(2))
 		{
 			u.cut_patch(i, rh.offset_abs(i, 0), rh.offset_abs(i, 1), rh.offset_abs(i, 2),
-						rh.size(i, 0), rh.size(i, 1), rh.size(i, 2));
+						rh.size(i, 0), rh.size(i, 1), rh.size(i, 2), benforce_coarse);
 		}
 	}
-
-	for (int i = rh.levelmax(); i > 0; --i)
-		mg_straight().restrict(*(u.get_grid(i)), *(u.get_grid(i - 1)));
 }

src/main.cc

@@ -456,7 +456,7 @@ int main(int argc, const char *argv[])
 	if (!cf.contains_key("setup", "kspace_TF"))
 		cf.insert_value("setup", "kspace_TF", "yes");
 
-	bool bspectral_sampling = cf.get_value<bool>("setup", "kspace_TF");
+	bool bspectral_sampling = true;//cf.get_value<bool>("setup", "kspace_TF");
 
 	if (bspectral_sampling)
 		music::ilog.Print("Using k-space sampled transfer functions...");

src/mesh.hh

@@ -8,8 +8,7 @@
 
 */
 
-#ifndef __MESH_HH
+#pragma once
-#define __MESH_HH
 
 #include <iostream>
 #include <iomanip>
@@ -27,7 +26,6 @@ using index_t = ptrdiff_t;
 using index3_t = std::array<index_t, 3>;
 using vec3_t = std::array<double, 3>;
 
-
 class refinement_mask
 {
 protected:
@@ -147,9 +145,11 @@ public:
 
 		m_pdata = new real_t[m_nx * m_ny * m_nz];
 
-		if (copy_over)
+		if (copy_over){
+			#pragma omp parallel for
 			for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
 				m_pdata[i] = m.m_pdata[i];
+		}
 	}
 
 	//! standard copy constructor
@@ -165,6 +165,7 @@ public:
 
 		m_pdata = new real_t[m_nx * m_ny * m_nz];
 
+		#pragma omp parallel for
 		for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
 			m_pdata[i] = m.m_pdata[i];
 	}
@@ -217,6 +218,7 @@ public:
 	//! set all the data to zero values
 	void zero(void)
 	{
+		#pragma omp parallel for
 		for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
 			m_pdata[i] = 0.0;
 	}
@@ -258,6 +260,7 @@ public:
 	//! direct multiplication of the whole data block with a number
 	Meshvar<real_t> &operator*=(real_t x)
 	{
+		#pragma omp parallel for
 		for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
 			m_pdata[i] *= x;
 		return *this;
@@ -266,6 +269,7 @@ public:
 	//! direct addition of a number to the whole data block
 	Meshvar<real_t> &operator+=(real_t x)
 	{
+		#pragma omp parallel for
 		for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
 			m_pdata[i] += x;
 		return *this;
@@ -274,6 +278,7 @@ public:
 	//! direct element-wise division of the whole data block by a number
 	Meshvar<real_t> &operator/=(real_t x)
 	{
+		#pragma omp parallel for
 		for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
 			m_pdata[i] /= x;
 		return *this;
@@ -282,6 +287,7 @@ public:
 	//! direct subtraction of a number from the whole data block
 	Meshvar<real_t> &operator-=(real_t x)
 	{
+		#pragma omp parallel for
 		for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
 			m_pdata[i] -= x;
 		return *this;
@@ -295,6 +301,7 @@ public:
 			music::elog.Print("Meshvar::operator*= : attempt to operate on incompatible data");
 			throw std::runtime_error("Meshvar::operator*= : attempt to operate on incompatible data");
 		}
+		#pragma omp parallel for
 		for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
 			m_pdata[i] *= v.m_pdata[i];
 
@@ -310,6 +317,7 @@ public:
 			throw std::runtime_error("Meshvar::operator/= : attempt to operate on incompatible data");
 		}
 
+		#pragma omp parallel for
 		for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
 			m_pdata[i] /= v.m_pdata[i];
 
@@ -324,6 +332,7 @@ public:
 			music::elog.Print("Meshvar::operator+= : attempt to operate on incompatible data");
 			throw std::runtime_error("Meshvar::operator+= : attempt to operate on incompatible data");
 		}
+		#pragma omp parallel for
 		for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
 			m_pdata[i] += v.m_pdata[i];
 
@@ -338,6 +347,7 @@ public:
 			music::elog.Print("Meshvar::operator-= : attempt to operate on incompatible data");
 			throw std::runtime_error("Meshvar::operator-= : attempt to operate on incompatible data");
 		}
+		#pragma omp parallel for
 		for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
 			m_pdata[i] -= v.m_pdata[i];
 
@@ -360,6 +370,7 @@ public:
 
 		m_pdata = new real_t[m_nx * m_ny * m_nz];
 
+		#pragma omp parallel for
 		for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
 			m_pdata[i] = m.m_pdata[i];
 
@@ -462,52 +473,13 @@ public:
 			m_pdata = new real_t[m_nx * m_ny * m_nz];
 		}
 
+		#pragma omp parallel for
 		for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
 			this->m_pdata[i] = m.m_pdata[i];
 
 		return *this;
 	}
 
-	//! sets the value of all ghost zones to zero
-	void zero_bnd(void)
-	{
-
-		int nx, ny, nz;
-		nx = this->size(0);
-		ny = this->size(1);
-		nz = this->size(2);
-
-		for (int j = -m_nbnd; j < ny + m_nbnd; ++j)
-			for (int k = -m_nbnd; k < nz + m_nbnd; ++k)
-			{
-				for (int i = -m_nbnd; i < 0; ++i)
-				{
-					(*this)(i, j, k) = 0.0;
-					(*this)(nx - 1 - i, j, k) = 0.0;
-				}
-			}
-
-		for (int i = -m_nbnd; i < nx + m_nbnd; ++i)
-			for (int k = -m_nbnd; k < nz + m_nbnd; ++k)
-			{
-				for (int j = -m_nbnd; j < 0; ++j)
-				{
-					(*this)(i, j, k) = 0.0;
-					(*this)(i, ny - j - 1, k) = 0.0;
-				}
-			}
-
-		for (int i = -m_nbnd; i < nx + m_nbnd; ++i)
-			for (int j = -m_nbnd; j < ny + m_nbnd; ++j)
-			{
-				for (int k = -m_nbnd; k < 0; ++k)
-				{
-					(*this)(i, j, k) = 0.0;
-					(*this)(i, j, nz - k - 1) = 0.0;
-				}
-			}
-	}
-
 	//! outputs the data, for debugging only, not practical for large datasets
 	void print(void) const
 	{
@@ -1121,8 +1093,9 @@ public:
 	 * @param nx new x-extent in fine grid cells
 	 * @param ny new y-extent in fine grid cells
 	 * @param nz new z-extent in fine grid cells
+	 * @param enforce_coarse_mean enforces the average of 8 cells on a fine level to equal the coarse 
 	 */
-	void cut_patch(unsigned ilevel, unsigned xoff, unsigned yoff, unsigned zoff, unsigned nx, unsigned ny, unsigned nz)
+	void cut_patch(unsigned ilevel, unsigned xoff, unsigned yoff, unsigned zoff, unsigned nx, unsigned ny, unsigned nz, bool enforce_coarse_mean)
 	{
 		unsigned dx, dy, dz, dxtop, dytop, dztop;
 
@@ -1139,10 +1112,19 @@ public:
 		MeshvarBnd<T> *mnew = new MeshvarBnd<T>(m_nbnd, nx, ny, nz, dxtop, dytop, dztop);
 
 		//... copy data
-		for (unsigned i = 0; i < nx; ++i)
+		[[maybe_unused]] double coarsesum = 0.0, finesum = 0.0;
-			for (unsigned j = 0; j < ny; ++j)
+    [[maybe_unused]] size_t coarsecount = 0, finecount = 0;
-				for (unsigned k = 0; k < nz; ++k)
+
-					(*mnew)(i, j, k) = (*m_pgrids[ilevel])(i + dx, j + dy, k + dz);
+    //... copy data
+		#pragma omp parallel for reduction(+:finesum,finecount) collapse(3)
+    for( unsigned i=0; i<nx; ++i )
+      for( unsigned j=0; j<ny; ++j )
+				for( unsigned k=0; k<nz; ++k )
+				{
+					(*mnew)(i,j,k) = (*m_pgrids[ilevel])(i+dx,j+dy,k+dz);
+					finesum += (*mnew)(i,j,k);
+					finecount++;
+				}
 
 		//... replace in hierarchy
 		delete m_pgrids[ilevel];
@@ -1160,6 +1142,37 @@ public:
 			m_pgrids[ilevel + 1]->offset(2) -= dz;
 		}
 
+		if( enforce_coarse_mean )
+		{
+			//... enforce top 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; ++i)
+					for (unsigned j = 0; j < ny; ++j)
+						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);
+			}
+		}
+
 		find_new_levelmin();
 	}
 
@@ -1207,14 +1220,14 @@ class refinement_hierarchy
 	std::vector<index3_t> len_;
 
 	unsigned
-			levelmin_,					//!< minimum grid level for Poisson solver
+			levelmin_,				//!< minimum grid level for Poisson solver
-			levelmax_,					//!< maximum grid level for all operations
+			levelmax_,				//!< maximum grid level for all operations
-			levelmin_tf_, 			//!< minimum grid level for density calculation
+			levelmin_tf_,			//!< minimum grid level for density calculation
-			padding_,						//!< padding in number of coarse cells between refinement levels
+			padding_,					//!< padding in number of coarse cells between refinement levels
-			blocking_factor_,		//!< blocking factor of grids, necessary fo BoxLib codes such as NyX
+			blocking_factor_, //!< blocking factor of grids, necessary fo BoxLib codes such as NyX
-			gridding_unit_;  	  //!< internal blocking factor of grids, necessary for Panphasia
+			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)
+	int margin_; //!< number of cells used for additional padding for convolutions with isolated boundaries (-1 = double padding)
 
 	config_file &cf_; //!< reference to config_file
 
@@ -1232,28 +1245,29 @@ class refinement_hierarchy
 	double rshift_[3];
 
 	//! calculates the greatest common divisor
-  int gcd(int a, int b) const {
+	int gcd(int a, int b) const
-    return b == 0 ? a : gcd(b, a % b);
+	{
-  }
+		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
+	// 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"
+	// 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 get_shift_unit(int base_unit, int levelmin) const
-    /*int Lp = 0;
+	{
-    while( base_unit * (1<<Lp) < 1<<(levelmin+1) ){
+		/*int Lp = 0;
-      ++Lp;
+		while( base_unit * (1<<Lp) < 1<<(levelmin+1) ){
-    }
+			++Lp;
-    int U = base_unit * (1<<Lp);
+		}
+		int U = base_unit * (1<<Lp);
 
-    return std::max<int>( 1, (1<<(levelmin+1)) / (2*gcd(U,1<<(levelmin+1) )) );*/
+		return std::max<int>( 1, (1<<(levelmin+1)) / (2*gcd(U,1<<(levelmin+1) )) );*/
 
-    int level_m = 0;
+		int level_m = 0;
-    while( base_unit * (1<<level_m) < (1<<levelmin) )
+		while (base_unit * (1 << level_m) < (1 << levelmin))
-      ++level_m;
+			++level_m;
 
-    return std::max<int>( 1, (1<<levelmin)/gcd(base_unit * (1<<level_m),(1<<levelmin)) );
+		return std::max<int>(1, (1 << levelmin) / gcd(base_unit * (1 << level_m), (1 << levelmin)));
-
+	}
-  }
 
 public:
 	//! copy constructor
@@ -1275,21 +1289,23 @@ public:
 		preserve_dims_ = cf_.get_value_safe<bool>("setup", "preserve_dims", false);
 		equal_extent_ = cf_.get_value_safe<bool>("setup", "force_equal_extent", false);
 		blocking_factor_ = cf.get_value_safe<unsigned>("setup", "blocking_factor", 0);
-		margin_          = cf.get_value_safe<int>("setup","convolution_margin",32);
+		margin_ = cf.get_value_safe<int>("setup", "convolution_margin", 32);
 
 		bool bnoshift = cf_.get_value_safe<bool>("setup", "no_shift", false);
 		bool force_shift = cf_.get_value_safe<bool>("setup", "force_shift", false);
 
 		gridding_unit_ = cf.get_value_safe<unsigned>("setup", "gridding_unit", 2);
 
-    if (gridding_unit_ != 2 && blocking_factor_==0) {
+		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_ ) {
+			blocking_factor_ = gridding_unit_; // THIS WILL LIKELY CAUSE PROBLEMS WITH NYX
-			music::elog.Print("incompatible gridding unit %d and blocking factor specified", gridding_unit_, blocking_factor_ );
+		}
+		else if (gridding_unit_ != 2 && blocking_factor_ != 0 && gridding_unit_ != blocking_factor_)
+		{
+			music::elog.Print("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_)
 		{
@@ -1312,8 +1328,8 @@ public:
 		// if not doing any refinement levels, set extent to full box
 		if (levelmin_ == levelmax_)
 		{
-			x0ref_ = { 0.0, 0.0, 0.0 };
+			x0ref_ = {0.0, 0.0, 0.0};
-			lxref_ = { 1.0, 1.0, 1.0 };
+			lxref_ = {1.0, 1.0, 1.0};
 		}
 
 		unsigned ncoarse = 1 << levelmin_;
@@ -1327,18 +1343,19 @@ public:
 
 		if ((levelmin_ != levelmax_) && (!bnoshift || force_shift))
 		{
-			int random_base_grid_unit = cf.get_value_safe<int>("random","base_unit",1);
+			int random_base_grid_unit = cf.get_value_safe<int>("random", "base_unit", 1);
-      int shift_unit = get_shift_unit( random_base_grid_unit, levelmin_ );
+			int shift_unit = get_shift_unit(random_base_grid_unit, levelmin_);
-      if( shift_unit != 1 ){
+			if (shift_unit != 1)
-        music::ilog.Print("volume can only be shifted by multiples of %d coarse cells.",shift_unit);
+			{
-      }
+				music::ilog.Print("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_[0] = (int)((0.5 - xc[0]) * (double)ncoarse / shift_unit + 0.5) * shift_unit; // ARJ(int)((0.5 - xc[0]) * ncoarse);
-      xshift_[2] = (int)((0.5-xc[2]) * (double)ncoarse / shift_unit + 0.5) * shift_unit;//ARJ(int)((0.5 - xc[2]) * 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_[1] = (int)((0.5 - xc[1]) * ncoarse);
-      // xshift_[2] = (int)((0.5 - xc[2]) * ncoarse);
+			// xshift_[2] = (int)((0.5 - xc[2]) * ncoarse);
 		}
 		else
 		{
@@ -1458,15 +1475,15 @@ public:
 		else
 		{
 			//... require alignment with coarser grid
-      music::ilog.Print("- Internal refinement bounding box: [%d,%d]x[%d,%d]x[%d,%d]", il, ir, jl, jr, kl, kr);
+			music::ilog.Print("- Internal refinement bounding box: [%d,%d]x[%d,%d]x[%d,%d]", il, ir, jl, jr, kl, kr);
 
-      il -= il % gridding_unit_;
+			il -= il % gridding_unit_;
-      jl -= jl % gridding_unit_;
+			jl -= jl % gridding_unit_;
-      kl -= kl % gridding_unit_;
+			kl -= kl % gridding_unit_;
 
-      ir = ((ir%gridding_unit_)!=0)? (ir/gridding_unit_ + 1)*gridding_unit_ : ir;
+			ir = ((ir % gridding_unit_) != 0) ? (ir / gridding_unit_ + 1) * gridding_unit_ : ir;
-      jr = ((jr%gridding_unit_)!=0)? (jr/gridding_unit_ + 1)*gridding_unit_ : jr;
+			jr = ((jr % gridding_unit_) != 0) ? (jr / gridding_unit_ + 1) * gridding_unit_ : jr;
-      kr = ((kr%gridding_unit_)!=0)? (kr/gridding_unit_ + 1)*gridding_unit_ : kr;
+			kr = ((kr % gridding_unit_) != 0) ? (kr / gridding_unit_ + 1) * gridding_unit_ : kr;
 		}
 
 		// if doing unigrid, set region to whole box
@@ -1500,8 +1517,7 @@ public:
 		{
 
 			absoffsets_[levelmax_] = {(il + nresmax) % nresmax, (jl + nresmax) % nresmax, (kl + nresmax) % nresmax};
-			len_[levelmax_] = { ir-il, jr-jl, kr-kl };
+			len_[levelmax_] = {ir - il, jr - jl, kr - kl};
-			
 
 			if (equal_extent_)
 			{
@@ -1512,16 +1528,16 @@ public:
 					throw std::runtime_error("Specified equal_extent=yes conflicting with ref_dims which are not equal.");
 				}
 				size_t ilevel = levelmax_;
-				index_t nmax = *std::max_element(len_[ilevel].begin(),len_[ilevel].end());
+				index_t nmax = *std::max_element(len_[ilevel].begin(), len_[ilevel].end());
 
-				index3_t dx = {0,0,0};
+				index3_t dx = {0, 0, 0};
-				for( int idim=0; idim<3; ++idim )
+				for (int idim = 0; idim < 3; ++idim)
 				{
 					dx[idim] = (index_t)((double)(nmax - len_[ilevel][idim]) * 0.5);
 					absoffsets_[ilevel][idim] -= dx[idim];
 					len_[ilevel][idim] = nmax;
 				}
-				
+
 				il = absoffsets_[ilevel][0];
 				jl = absoffsets_[ilevel][1];
 				kl = absoffsets_[ilevel][2];
@@ -1575,12 +1591,12 @@ public:
 			{
 				//... require alignment with coarser grid
 				il -= il % gridding_unit_;
-        jl -= jl % gridding_unit_;
+				jl -= jl % gridding_unit_;
-        kl -= kl % gridding_unit_;
+				kl -= kl % gridding_unit_;
 
-        ir = ((ir%gridding_unit_)!=0)? (ir/gridding_unit_ + 1)*gridding_unit_ : ir;
+				ir = ((ir % gridding_unit_) != 0) ? (ir / gridding_unit_ + 1) * gridding_unit_ : ir;
-        jr = ((jr%gridding_unit_)!=0)? (jr/gridding_unit_ + 1)*gridding_unit_ : jr;
+				jr = ((jr % gridding_unit_) != 0) ? (jr / gridding_unit_ + 1) * gridding_unit_ : jr;
-        kr = ((kr%gridding_unit_)!=0)? (kr/gridding_unit_ + 1)*gridding_unit_ : kr;
+				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)
@@ -1588,19 +1604,19 @@ public:
 				music::elog.Print("Internal refinement bounding box error: [%d,%d]x[%d,%d]x[%d,%d], level=%d", il, ir, jl, jr, kl, kr, ilevel);
 				throw std::runtime_error("refinement_hierarchy: Internal refinement bounding box error 2");
 			}
-			absoffsets_[ilevel] = { il, jl, kl };
+			absoffsets_[ilevel] = {il, jl, kl};
-			len_[ilevel] = { ir-il, jr-jl, kr-kl };
+			len_[ilevel] = {ir - il, jr - jl, kr - kl};
 
 			if (blocking_factor_)
 			{
-				for( int idim=0; idim<3; ++idim )
+				for (int idim = 0; idim < 3; ++idim)
-					len_[ilevel][idim] += len_[ilevel][idim] % blocking_factor_; 
+					len_[ilevel][idim] += len_[ilevel][idim] % blocking_factor_;
 			}
 
 			if (equal_extent_)
 			{
-				index_t nmax = *std::max_element(len_[ilevel].begin(),len_[ilevel].end());
+				index_t nmax = *std::max_element(len_[ilevel].begin(), len_[ilevel].end());
-				for( int idim=0; idim<3; ++idim )
+				for (int idim = 0; idim < 3; ++idim)
 				{
 					index_t dx = (int)((double)(nmax - len_[ilevel][idim]) * 0.5);
 					absoffsets_[ilevel][idim] -= dx;
@@ -1618,17 +1634,17 @@ public:
 
 		//... determine relative offsets between grids
 		for (unsigned ilevel = levelmax_; ilevel > levelmin_; --ilevel)
-			for(int idim=0; idim<3; ++idim )
+			for (int idim = 0; idim < 3; ++idim)
 			{
-				offsets_[ilevel][idim] = (absoffsets_[ilevel][idim]/2 - absoffsets_[ilevel-1][idim]);
+				offsets_[ilevel][idim] = (absoffsets_[ilevel][idim] / 2 - absoffsets_[ilevel - 1][idim]);
 			}
 
 		//... do a forward sweep to ensure that absolute offsets are also correct now
 		for (unsigned ilevel = levelmin_ + 1; ilevel <= levelmax_; ++ilevel)
-			for(int idim=0; idim<3; ++idim )
+			for (int idim = 0; idim < 3; ++idim)
-				{
+			{
-					absoffsets_[ilevel][idim] = 2 * absoffsets_[ilevel-1][idim] + 2 * offsets_[ilevel][idim];
+				absoffsets_[ilevel][idim] = 2 * absoffsets_[ilevel - 1][idim] + 2 * offsets_[ilevel][idim];
-				}
+			}
 
 		for (unsigned ilevel = levelmin_ + 1; ilevel <= levelmax_; ++ilevel)
 		{
@@ -1751,7 +1767,7 @@ public:
 		{
 			unsigned n = 1 << i;
 
-			if ( absoffsets_[i] == index3_t({0,0,0}) && len_[i]== index3_t({n,n,n}) )
+			if (absoffsets_[i] == index3_t({0, 0, 0}) && len_[i] == index3_t({n, n, n}))
 			{
 				levelmin_ = i;
 			}
@@ -1816,7 +1832,7 @@ public:
 
 		if (xshift_[0] != 0 || xshift_[1] != 0 || xshift_[2] != 0)
 			music::ilog << " - Domain will be shifted by (" << xshift_[0] << ", " << xshift_[1] << ", " << xshift_[2] << ")\n"
-								<< std::endl;
+									<< std::endl;
 
 		music::ilog << " - Grid structure:\n";
 
@@ -1844,5 +1860,3 @@ public:
 typedef GridHierarchy<real_t> grid_hierarchy;
 typedef MeshvarBnd<real_t> meshvar_bnd;
 typedef Meshvar<real_t> meshvar;
-
-#endif

src/mg_solver.hh

@@ -192,7 +192,7 @@ void solver<S, I, O>::GaussSeidel(real_t h, MeshvarBnd<real_t> *u, const Meshvar
 			h2 = h * h;
 
 	for (int color = 0; color < 2; ++color)
-#pragma omp parallel for
+		#pragma omp parallel for collapse(3)
 		for (int ix = 0; ix < nx; ++ix)
 			for (int iy = 0; iy < ny; ++iy)
 				for (int iz = 0; iz < nz; ++iz)
@@ -279,7 +279,7 @@ void solver<S, I, O>::twoGrid(unsigned ilevel)
 			ozp = uf->offset(2);
 
 	meshvar_bnd tLu(*uc, false);
-#pragma omp parallel for
+	#pragma omp parallel for
 	for (int ix = 0; ix < nx / 2; ++ix)
 	{
 		int iix = 2 * ix;
@@ -297,7 +297,7 @@ void solver<S, I, O>::twoGrid(unsigned ilevel)
 	oj = ff->offset(1);
 	ok = ff->offset(2);
 
-#pragma omp parallel for
+#pragma omp parallel for collapse(3)
 	for (int ix = oi; ix < oi + (int)ff->size(0) / 2; ++ix)
 		for (int iy = oj; iy < oj + (int)ff->size(1) / 2; ++iy)
 			for (int iz = ok; iz < ok + (int)ff->size(2) / 2; ++iz)
@@ -319,7 +319,7 @@ void solver<S, I, O>::twoGrid(unsigned ilevel)
 	meshvar_bnd cc(*uc, false);
 
 //... compute correction on coarse grid
-#pragma omp parallel for
+#pragma omp parallel for collapse(3)
 	for (int ix = 0; ix < (int)cc.size(0); ++ix)
 		for (int iy = 0; iy < (int)cc.size(1); ++iy)
 			for (int iz = 0; iz < (int)cc.size(2); ++iz)
@@ -372,8 +372,7 @@ double solver<S, I, O>::compute_error(const MeshvarBnd<real_t> &u, const Meshvar
 
 	double h = 1.0 / (1ul << ilevel), h2 = h * h;
 
-#pragma omp parallel for reduction(+ \
+#pragma omp parallel for reduction(+ : err, count) collapse(3)
-																	 : err, count)
 	for (int ix = 0; ix < nx; ++ix)
 		for (int iy = 0; iy < ny; ++iy)
 			for (int iz = 0; iz < nz; ++iz)
@@ -410,8 +409,7 @@ double solver<S, I, O>::compute_error(const GridHierarchy<real_t> &uh, const Gri
 
 		double h = 1.0 / (1ul << ilevel), h2 = h * h;
 
-#pragma omp parallel for reduction(+ \
+		#pragma omp parallel for reduction(+ : err, count, mean_res) collapse(3)
-																	 : err, count)
 		for (int ix = 0; ix < nx; ++ix)
 			for (int iy = 0; iy < ny; ++iy)
 				for (int iz = 0; iz < nz; ++iz)
@@ -461,8 +459,7 @@ double solver<S, I, O>::compute_RMS_resid(const GridHierarchy<real_t> &uh, const
 		double sum = 0.0, sumd2 = 0.0;
 		size_t count = 0;
 
-#pragma omp parallel for reduction(+ \
+		#pragma omp parallel for reduction(+ : sum, sumd2, count)
-																	 : sum, sumd2, count)
 		for (int ix = 0; ix < nx; ++ix)
 			for (int iy = 0; iy < ny; ++iy)
 				for (int iz = 0; iz < nz; ++iz)

src/plugins/random_music.cc

@@ -531,22 +531,16 @@ void RNG_music::fill_grid(int ilevel, DensityGrid<real_t> &A)
 
       if (nx == (int)A.size(0) + 2 * margin[0] && ny == (int)A.size(1) + 2 * margin[1] && nz == (int)A.size(2) + 2 * margin[2])
       {
-        int ox = margin[0], oy = margin[1], oz = margin[2];
         std::vector<real_t> slice(ny * nz, 0.0);
 
         for (int i = 0; i < nx; ++i)
         {
           ifs.read(reinterpret_cast<char *>(&slice[0]), ny * nz * sizeof(real_t));
 
-          if (i < ox)
-            continue;
-          if (i >= 3 * ox)
-            break;
-
           #pragma omp parallel for
-          for (int j = oy; j < 3 * oy; ++j)
+          for (int j = 0; j < ny; ++j)
-            for (int k = oz; k < 3 * oz; ++k)
+            for (int k = 0; k < nz; ++k)
-              A(i - ox, j - oy, k - oz) = slice[j * nz + k];
+              A(i, j, k) = slice[j * nz + k];
         }
 
         ifs.close();
@@ -577,7 +571,7 @@ void RNG_music::fill_grid(int ilevel, DensityGrid<real_t> &A)
   }
   else
   {
-    music::ulog.Print("Copying white noise from memory cache...");
+    music::ilog.Print("Copying white noise from memory cache...");
 
     if (mem_cache_[ilevel - levelmin_] == NULL)
       music::elog.Print("Tried to access mem-cached random numbers for level %d. But these are not available!\n", ilevel);

src/plugins/random_music_wnoise_generator.cc

@@ -627,16 +627,18 @@ music_wnoise_generator<T>::music_wnoise_generator(music_wnoise_generator<T> &rc,
 
         val *= val_phas * sqrt8;
 
-        if (i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2)
+        // if (x0_ == NULL || lx_ == NULL){
-        {
+          if (i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2)
-          RE(cfine[qf]) = val.real();
+          {
-          IM(cfine[qf]) = val.imag();
+            RE(cfine[qf]) = val.real();
-        }
+            IM(cfine[qf]) = val.imag();
-        else
+          }
-        {
+          else
-          // RE(cfine[qf]) = val.real();
+          {
-          // IM(cfine[qf]) = 0.0;
+            RE(cfine[qf]) = val.real();
-        }
+            IM(cfine[qf]) = 0.0;
+          }
+        // }
       }
 
     delete[] rcoarse;

src/plugins/random_panphasia.cc

@@ -223,10 +223,9 @@ public:
       incongruent_fields_ = true;
       ngrid_ = 2 * ratio * pdescriptor_->i_base;
       grid_rescale_fac_ = (double)ngrid_ / (1 << levelmin_);
-      music::ilog.Print("PANPHASIA: will use a higher resolution:\n"
+      music::ilog << "PANPHASIA: will use a higher resolution:" << std::endl
-              "     (%d -> %d) * 2**ref compatible with PANPHASIA\n"
+             << "     (" << grid_m_ << " -> " << grid_p_ 
-              "     will Fourier interpolate after.",
+             << ") * 2**ref compatible with PANPHASIA (will Fourier interpolate after)" << std::endl;
-              grid_m_, grid_p_);
     }
   }
 
@@ -790,7 +789,9 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
   delete[] pr3;
   delete[] pr4;
 
-  music::ulog.Print("Copying random field data %d,%d,%d -> %d,%d,%d", nxremap[0], nxremap[1], nxremap[2], nx[0], nx[1], nx[2]);
+  music::ilog.Print("Copying random field data %d,%d,%d -> %d,%d,%d", nxremap[0], nxremap[1], nxremap[2], nx[0], nx[1], nx[2]);
+  music::ilog.Print("iexpand_levt = %d,%d,%d", iexpand_left[0], iexpand_left[1], iexpand_left[2]);
+  
 
   //    n = 1<<level;
   //    ng = n;
@@ -829,7 +830,8 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
   music::ulog.Print("done with PANPHASIA for level %d:\n       mean=%g, var=%g", level, sum, sum2);
   music::ulog.Print("Copying into R array: nx[0],nx[1],nx[2] %d %d %d \n", nx[0], nx[1], nx[2]);
 
-  music::ilog.Print("PANPHASIA level %d mean and variance are\n       <p> = %g | var(p) = %g", level, sum, sum2);
+  music::ilog << "PANPHASIA level " << level << " mean and variance are" << std::endl
+              << "       <p> = "<< sum << " | var(p) = " << sum2 << std::endl;
 }
 
 namespace