From 3ad0ac6a108321dbf219d4075089aa6f11a195f4 Mon Sep 17 00:00:00 2001
From: Oliver Hahn <oliverjhahn@gmail.com>
Date: Tue, 16 Nov 2010 20:56:08 -0800
Subject: [PATCH] Fixed a memory leak. Slightly better memory conservation.
 Fixed a bug introduced in the previous revision that led to code termination.

---
 convolution_kernel.cc |  2 +-
 main.cc               | 71 ++++++++++++++++++++++++++++++++++---------
 mesh.hh               | 20 +++++++-----
 mg_solver.hh          | 30 ++++++++++++++++--
 poisson.cc            | 19 +++++++-----
 5 files changed, 111 insertions(+), 31 deletions(-)

diff --git a/convolution_kernel.cc b/convolution_kernel.cc
index f522b5c..b04562e 100644
--- a/convolution_kernel.cc
+++ b/convolution_kernel.cc
@@ -70,7 +70,7 @@ namespace convolution{
 		//..... need a phase shift for baryons for SPH
 		bool do_SPH = pk->pcf_->getValueSafe<bool>("setup","do_SPH",false) & (pk->type_==baryon | pk->type_==vbaryon);
 		double dsph = 0.0;
-		double boxlength = pk->pcf_->getValue<int>("setup","boxlength");
+		double boxlength = pk->pcf_->getValue<double>("setup","boxlength");
 		if( do_SPH )
 		{	
 			int lmax = pk->pcf_->getValue<int>("setup","levelmax");
diff --git a/main.cc b/main.cc
index fdbefac..b92e1e0 100644
--- a/main.cc
+++ b/main.cc
@@ -416,12 +416,28 @@ int main (int argc, const char * argv[])
 	//------------------------------------------------------------------------------
 	//... initialize the random numbers
 	//------------------------------------------------------------------------------
+	std::cout << "=============================================================\n";
+	std::cout << "   GENERATING WHITE NOISE\n";
+	std::cout << "-------------------------------------------------------------\n";
+	LOGUSER("Computing white noise...");
 	rand_gen rand( cf, rh_TF, the_transfer_function_plugin );
 	
 	//------------------------------------------------------------------------------
 	//... initialize the Poisson solver
 	//------------------------------------------------------------------------------
+	bool bdefd = cf.getValueSafe<bool> ( "poisson" , "fft_fine", true );
+	bool bsph  = cf.getValueSafe<bool>("setup","do_SPH",false) && do_baryons;
+	
+	
 	bool kspace	= cf.getValueSafe<bool>( "poisson", "kspace", false );
+	
+	//... if in unigrid mode, use k-space instead of hybrid
+	if(bdefd&lbase==lmax)
+	{
+		kspace=true;
+		bdefd=false;
+	}
+	
 	std::string poisson_solver_name;
 	if( kspace )
 		poisson_solver_name = std::string("fft_poisson");
@@ -429,15 +445,12 @@ int main (int argc, const char * argv[])
 		poisson_solver_name = std::string("mg_poisson");
 	
 	unsigned grad_order = cf.getValueSafe<unsigned> ( "poisson" , "grad_order", 4 );
-	bool bdefd = cf.getValueSafe<bool> ( "poisson" , "fft_fine", true );
-	bool bsph  = cf.getValueSafe<bool>("setup","do_SPH",false) && do_baryons;
 	
-	//... if in unigrid mode, use k-space instead
-	//if(bdefd&lbase==lmax)
-	//kspace=true;
+	
+	
 	
 	//... switch off if using kspace anyway
-	bdefd &= !kspace;
+	//bdefd &= !kspace;
 	
 	poisson_plugin_creator *the_poisson_plugin_creator = get_poisson_plugin_map()[ poisson_solver_name ];
 	poisson_plugin *the_poisson_solver = the_poisson_plugin_creator->create( cf );
@@ -459,7 +472,7 @@ int main (int argc, const char * argv[])
 			std::cout << "-------------------------------------------------------------\n";
 			LOGUSER("Computing dark matter displacements...");
 			
-			grid_hierarchy f( nbnd ), u(nbnd);
+			grid_hierarchy f( nbnd );//, u(nbnd);
 			tf_type my_tf_type = cdm;
 			if( !do_baryons || !the_transfer_function_plugin->tf_is_distinct() )
 				my_tf_type = total;
@@ -469,11 +482,10 @@ int main (int argc, const char * argv[])
 			coarsen_density(rh_Poisson, f);
 			normalize_density(f);
 			
-			u = f;	u.zero();
-			
 			the_output_plugin->write_dm_mass(f);
 			the_output_plugin->write_dm_density(f);
 			
+			grid_hierarchy u( f );	u.zero();
 			err = the_poisson_solver->solve(f, u);
 			
 			if(!bdefd)
@@ -504,8 +516,11 @@ int main (int argc, const char * argv[])
 
 					the_output_plugin->write_dm_position(icoord, data_forIO );
 				}
+				u.deallocate();
+				data_forIO.deallocate();
 			}
 			
+			
 			//------------------------------------------------------------------------------
 			//... gas density
 			//------------------------------------------------------------------------------
@@ -526,6 +541,9 @@ int main (int argc, const char * argv[])
 					u = f;	u.zero();
 					err = the_poisson_solver->solve(f, u);
 					
+					if(!bdefd)
+						f.deallocate();
+					
 					grid_hierarchy data_forIO(u);
 					for( int icoord = 0; icoord < 3; ++icoord )
 					{
@@ -542,18 +560,25 @@ int main (int argc, const char * argv[])
 							//... displacement
 							the_poisson_solver->gradient(icoord, u, data_forIO );
 						
+						
 						the_output_plugin->write_gas_position(icoord, data_forIO );
+						
 					}
+					u.deallocate();
+					data_forIO.deallocate();
+					f.deallocate();
 				}
 				else if( do_LLA )
 				{
 					u = f;	u.zero();
 					err = the_poisson_solver->solve(f, u);
 					compute_LLA_density( u, f,grad_order );
+					u.deallocate();
 					normalize_density(f);
 				}
 				
 				the_output_plugin->write_gas_density(f);
+				f.deallocate();
 			}
 			
 			
@@ -596,6 +621,8 @@ int main (int argc, const char * argv[])
 					else 
 						the_poisson_solver->gradient(icoord, u, data_forIO );
 					
+					
+					
 					//... multiply to get velocity
 					data_forIO *= cosmo.vfact;
 					
@@ -606,8 +633,13 @@ int main (int argc, const char * argv[])
 
 					if( do_baryons )
 						the_output_plugin->write_gas_velocity(icoord, data_forIO);
+					
+				
 				}
 				
+				u.deallocate();
+				data_forIO.deallocate();
+				
 			}
 			else
 			{
@@ -655,7 +687,9 @@ int main (int argc, const char * argv[])
 					
 					the_output_plugin->write_dm_velocity(icoord, data_forIO);
 				}
+				u.deallocate();
 				data_forIO.deallocate();
+				f.deallocate();
 				
 				
 				std::cout << "=============================================================\n";
@@ -698,6 +732,9 @@ int main (int argc, const char * argv[])
 					
 					the_output_plugin->write_gas_velocity(icoord, data_forIO);
 				}
+				u.deallocate();
+				f.deallocate();
+				data_forIO.deallocate();
 			}
 		/*********************************************************************************************/
 		/*********************************************************************************************/
@@ -794,6 +831,8 @@ int main (int argc, const char * argv[])
 					the_output_plugin->write_gas_velocity(icoord, data_forIO);				
 			}
 			data_forIO.deallocate();
+			if( !dm_only )
+				u1.deallocate();
 			
 			
 			if( do_baryons && the_transfer_function_plugin->tf_has_velocities() )
@@ -864,6 +903,7 @@ int main (int argc, const char * argv[])
 					the_output_plugin->write_gas_velocity(icoord, data_forIO);				
 				}
 				data_forIO.deallocate();
+				u1.deallocate();
 			}
 			
 			
@@ -952,6 +992,9 @@ int main (int argc, const char * argv[])
 				the_output_plugin->write_dm_position(icoord, data_forIO );	
 			}
 			
+			data_forIO.deallocate();
+			u1.deallocate();
+			
 
 			if( do_baryons && !bsph )
 			{	
@@ -1054,7 +1097,12 @@ int main (int argc, const char * argv[])
 
 		}
 		
+		//------------------------------------------------------------------------------
+		//... finish output
+		//------------------------------------------------------------------------------
 		
+		the_output_plugin->finalize();
+		delete the_output_plugin;
 		
 	}catch(std::runtime_error& excp){
 		LOGERR("Fatal error occured. Code will exit:");
@@ -1066,12 +1114,7 @@ int main (int argc, const char * argv[])
 
 	std::cout << "=============================================================\n";
 	
-	//------------------------------------------------------------------------------
-	//... finish output
-	//------------------------------------------------------------------------------
 	
-	the_output_plugin->finalize();
-	delete the_output_plugin;
 
 	if( !bfatal )
 	{	
diff --git a/mesh.hh b/mesh.hh
index e39f0c0..8ee8ccc 100644
--- a/mesh.hh
+++ b/mesh.hh
@@ -337,14 +337,17 @@ public:
 	//! assignment operator for rectangular meshes with ghost zones
 	MeshvarBnd<real_t>& operator=( const MeshvarBnd<real_t>& m )
 	{
-		this->m_nx = m.m_nx;
-		this->m_ny = m.m_ny;
-		this->m_nz = m.m_nz;
-		
-		if( m_pdata != NULL )
-			delete[] m_pdata;
+		if( this->m_nx != m.m_nx || this->m_ny != m.m_ny || this->m_nz != m.m_nz )
+		{
+			this->m_nx = m.m_nx;
+			this->m_ny = m.m_ny;
+			this->m_nz = m.m_nz;
 			
-		m_pdata = new real_t[m_nx*m_ny*m_nz];
+			if( m_pdata != NULL )
+				delete[] m_pdata;
+			
+			m_pdata = new real_t[m_nx*m_ny*m_nz];			
+		}
 		
 		for( size_t i=0; i<m_nx*m_ny*m_nz; ++i )
 			this->m_pdata[i] = m.m_pdata[i];
@@ -665,6 +668,9 @@ public:
 	void create_base_hierarchy( unsigned lmax )
 	{
 		unsigned n=1;
+		
+		this->deallocate();
+		
 		m_pgrids.clear();
 		
 		m_xoffabs.clear();
diff --git a/mg_solver.hh b/mg_solver.hh
index ca3d167..f7ece6d 100644
--- a/mg_solver.hh
+++ b/mg_solver.hh
@@ -296,7 +296,7 @@ void solver<S,I,O,T>::twoGrid( unsigned ilevel )
 		
 	
 
-	meshvar_bnd Lu(*uf,false);
+	/*meshvar_bnd Lu(*uf,false);
 	Lu.zero();
 
 	#pragma omp parallel for
@@ -310,7 +310,33 @@ void solver<S,I,O,T>::twoGrid( unsigned ilevel )
 	
 	//... restrict Lu
 	m_gridop.restrict( Lu, tLu );
-	Lu.deallocate();
+	Lu.deallocate();*/
+	
+	int 
+	oxp = uf->offset(0),
+	oyp = uf->offset(1),
+	ozp = uf->offset(2);
+	
+	meshvar_bnd tLu(*uc,false);
+	#pragma omp parallel for
+	for( int ix=0; ix<nx/2; ++ix )
+	{	
+		int iix=2*ix;
+		for( int iy=0,iiy=0; iy<ny/2; ++iy,iiy+=2 )
+		
+		
+			for( int iz=0,iiz=0; iz<nz/2; ++iz,iiz+=2 )
+				tLu(ix+oxp,iy+oyp,iz+ozp) = 0.125 * (
+							 m_scheme.apply( (*uf), iix, iiy, iiz )
+							+m_scheme.apply( (*uf), iix, iiy, iiz+1 )
+							+m_scheme.apply( (*uf), iix, iiy+1, iiz )
+							+m_scheme.apply( (*uf), iix, iiy+1, iiz+1 )
+							+m_scheme.apply( (*uf), iix+1, iiy, iiz )
+							+m_scheme.apply( (*uf), iix+1, iiy, iiz+1 )
+							+m_scheme.apply( (*uf), iix+1, iiy+1, iiz )
+							+m_scheme.apply( (*uf), iix+1, iiy+1, iiz+1 )
+						)/h2;
+	}
 	
 	//... restrict source term
 	m_gridop.restrict( *ff, *fc );
diff --git a/poisson.cc b/poisson.cc
index effb500..9b575ac 100644
--- a/poisson.cc
+++ b/poisson.cc
@@ -500,11 +500,12 @@ double fft_poisson_plugin::solve( grid_hierarchy& f, grid_hierarchy& u )
 	fftw_real *data = new fftw_real[nx*ny*nzp];
 	fftw_complex *cdata = reinterpret_cast<fftw_complex*> (data);
 	
+	#pragma omp parallel for
 	for( int i=0; i<nx; ++i )
 		for( int j=0; j<ny; ++j )	
 			for( int k=0; k<nz; ++k )
 			{
-				unsigned idx = (i*ny+j)*nzp+k;
+				size_t idx = ((size_t)i*ny+(size_t)j)*nzp+(size_t)k;
 				data[idx] = (*f.get_grid(f.levelmax()))(i,j,k);
 			}
 	
@@ -536,6 +537,7 @@ double fft_poisson_plugin::solve( grid_hierarchy& f, grid_hierarchy& u )
 	double kfac = 2.0*M_PI;
 	double fac = -1.0/(nx*ny*nz);
 	
+	#pragma omp parallel for
 	for( int i=0; i<nx; ++i )
 		for( int j=0; j<ny; ++j )	
 			for( int k=0; k<nz/2+1; ++k )
@@ -548,7 +550,7 @@ double fft_poisson_plugin::solve( grid_hierarchy& f, grid_hierarchy& u )
 				
 				double kk2 = kfac*kfac*(ki*ki+kj*kj+kk*kk);
 				
-				unsigned idx = (i*ny+j)*nzp/2+k;
+				size_t idx = ((size_t)i*ny+(size_t)j)*nzp/2+(size_t)k;
 #ifdef FFTW3
 				cdata[idx][0] *= -1.0/kk2*fac;
 				cdata[idx][1] *= -1.0/kk2*fac;
@@ -588,11 +590,12 @@ double fft_poisson_plugin::solve( grid_hierarchy& f, grid_hierarchy& u )
 
 	
 	//... copy data ..........................................
+	#pragma omp parallel for
 	for( int i=0; i<nx; ++i )
 		for( int j=0; j<ny; ++j )	
 			for( int k=0; k<nz; ++k )
 			{
-				unsigned idx = (i*ny+j)*nzp+k;
+				size_t idx = ((size_t)i*ny+(size_t)j)*nzp+(size_t)k;
 				(*u.get_grid(u.levelmax()))(i,j,k) = data[idx];
 			}
 	
@@ -622,11 +625,12 @@ double fft_poisson_plugin::gradient( int dir, grid_hierarchy& u, grid_hierarchy&
 	fftw_real *data = new fftw_real[nx*ny*nzp];
 	fftw_complex *cdata = reinterpret_cast<fftw_complex*> (data);
 	
+	#pragma omp parallel for
 	for( int i=0; i<nx; ++i )
 		for( int j=0; j<ny; ++j )	
 			for( int k=0; k<nz; ++k )
 			{
-				unsigned idx = (i*ny+j)*nzp+k;
+				size_t idx = ((size_t)i*ny+(size_t)j)*nzp+(size_t)k;
 				data[idx] = (*u.get_grid(u.levelmax()))(i,j,k);
 			}
 	
@@ -658,11 +662,12 @@ double fft_poisson_plugin::gradient( int dir, grid_hierarchy& u, grid_hierarchy&
 	double fac = -1.0/(nx*ny*nz);
 	double kfac = 2.0*M_PI;
 	
+	#pragma omp parallel for
 	for( int i=0; i<nx; ++i )
 		for( int j=0; j<ny; ++j )	
 			for( int k=0; k<nz/2+1; ++k )
 			{
-				unsigned idx = (i*ny+j)*nzp/2+k;
+				size_t idx = (i*(size_t)ny+j)*(size_t)nzp/2+(size_t)k;
 				int ii = i; if(ii>nx/2) ii-=nx;
 				int jj = j; if(jj>ny/2) jj-=ny;
 				double ki = (double)ii;
@@ -727,7 +732,7 @@ double fft_poisson_plugin::gradient( int dir, grid_hierarchy& u, grid_hierarchy&
 		for( int j=0; j<ny; ++j )	
 			for( int k=0; k<nz; ++k )
 			{
-				unsigned idx = (i*ny+j)*nzp+k;
+				size_t idx = ((size_t)i*ny+(size_t)j)*nzp+(size_t)k;
 				(*Du.get_grid(u.levelmax()))(i,j,k) = data[idx];
 				if(fabs(data[idx])>dmax)
 					dmax = fabs(data[idx]);
@@ -1023,7 +1028,7 @@ void poisson_hybrid( MeshvarBnd<double>& f, int idir, int order, bool periodic )
 	
 	size_t N = (size_t)nxp*(size_t)nyp*2*((size_t)nzp/2+1);
 	
-	//#pragma omp parallel for
+	#pragma omp parallel for
 	for( size_t i=0; i<N; ++i )
 		data[i]=0.0;