mirror of
https://github.com/cosmo-sims/MUSIC.git
synced 2024-09-19 17:03:46 +02:00
Fixed a memory leak.
Slightly better memory conservation. Fixed a bug introduced in the previous revision that led to code termination.
This commit is contained in:
Oliver Hahn
parent
8d9068852e
commit
3ad0ac6a10
5 changed files with 111 additions and 31 deletions
|
|
@ -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");
|
||||
|
|
|
|||
71
main.cc
71
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 )
|
||||
{
|
||||
|
|
|
|||
20
mesh.hh
20
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();
|
||||
|
|
|
|||
30
mg_solver.hh
30
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 );
|
||||
|
|
|
|||
19
poisson.cc
19
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;
|
||||
|
||||
|
|
|
|||
Loading…
Reference in a new issue