diff --git a/densities.cc b/densities.cc index bba3c09..31ec217 100644 --- a/densities.cc +++ b/densities.cc @@ -4,7 +4,7 @@ a code to generate multi-scale initial conditions for cosmological simulations - Copyright (C) 2010-13 Oliver Hahn + Copyright (C) 2010 Oliver Hahn */ @@ -27,9 +27,9 @@ void fft_interpolate( m1& V, m2& v, bool from_basegrid=false ) if( !from_basegrid ) { - oxf += nxF/4; - oyf += nyF/4; - ozf += nzF/4; + oxf += nxF/4; + oyf += nyF/4; + ozf += nzF/4; } LOGUSER("FFT interpolate: offset=%d,%d,%d size=%d,%d,%d",oxf,oyf,ozf,nxf,nyf,nzf); @@ -58,16 +58,15 @@ void fft_interpolate( m1& V, m2& v, bool from_basegrid=false ) size_t q = ((size_t)ii*nyc+(size_t)jj)*nzcp+(size_t)kk; rcoarse[q] = V( oxf+i, oyf+j, ozf+k ); } - + #pragma omp parallel for for( int i=0; i<(int)nxf; ++i ) for( int j=0; j<(int)nyf; ++j ) for( int k=0; k<(int)nzf; ++k ) { size_t q = ((size_t)i*nyf+(size_t)j)*nzfp+(size_t)k; - rfine[q] = v(i,j,k); + rfine[q] = v(i,j,k); } - #ifdef FFTW3 #ifdef SINGLE_PRECISION @@ -244,7 +243,6 @@ void fft_interpolate( m1& V, m2& v, bool from_basegrid=false ) } delete[] rfine; - } @@ -274,7 +272,7 @@ void GenerateDensityUnigrid( config_file& cf, transfer_function *ptf, tf_type ty if( kspace ) { - std::cout << " - Using new k-space transfer function kernel.\n"; + std::cout << " - Using k-space transfer function kernel.\n"; LOGUSER("Using k-space transfer function kernel."); #ifdef SINGLE_PRECISION @@ -334,84 +332,63 @@ void GenerateDensityUnigrid( config_file& cf, transfer_function *ptf, tf_type ty /*******************************************************************************************/ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type type, - refinement_hierarchy& refh, rand_gen& rand, grid_hierarchy& delta, bool smooth, bool shift ) + refinement_hierarchy& refh, rand_gen& rand, + grid_hierarchy& delta, bool smooth, bool shift ) { - unsigned levelmin,levelmax,levelminPoisson; - std::vector rngseeds; - std::vector rngfnames; - bool kspaceTF; - - double tstart, tend; - + unsigned levelmin, levelmax, levelminPoisson; + std::vector rngseeds; + std::vector rngfnames; + bool kspaceTF; + + double tstart, tend; + #ifndef SINGLETHREAD_FFTW - tstart = omp_get_wtime(); + tstart = omp_get_wtime(); #else - tstart = (double)clock() / CLOCKS_PER_SEC; + tstart = (double)clock() / CLOCKS_PER_SEC; #endif + + levelminPoisson = cf.getValue("setup","levelmin"); + levelmin = cf.getValueSafe("setup","levelmin_TF",levelminPoisson); + levelmax = cf.getValue("setup","levelmax"); + kspaceTF = cf.getValueSafe("setup", "kspace_TF", false); + + unsigned nbase = 1<("setup","levelmin"); - levelmin = cf.getValueSafe("setup","levelmin_TF",levelminPoisson); - levelmax = cf.getValue("setup","levelmax"); - kspaceTF = cf.getValueSafe("setup", "kspace_TF", false); + convolution::kernel_creator *the_kernel_creator; - - unsigned nbase = 1<create( cf, ptf, refh, type ); + convolution::kernel *the_tf_kernel = the_kernel_creator->create( cf, ptf, refh, type ); - - /***** PERFORM CONVOLUTIONS *****/ - - + /***** PERFORM CONVOLUTIONS *****/ if( kspaceTF ){ //... create and initialize density grids with white noise DensityGrid *top(NULL); PaddedDensitySubGrid *coarse(NULL), *fine(NULL); int nlevels = (int)levelmax-(int)levelmin+1; - + // do coarse level top = new DensityGrid( nbase, nbase, nbase ); LOGINFO("Performing noise convolution on level %3d",levelmin); @@ -424,34 +401,44 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type for( int i=1; i(refh.offset(levelmin+i,0), refh.offset(levelmin+i,1), refh.offset(levelmin+i,2), - refh.size(levelmin+i,0), refh.size(levelmin+i,1), refh.size(levelmin+i,2) ); + LOGUSER(" offset=(%5d,%5d,%5d)",refh.offset(levelmin+i,0), + refh.offset(levelmin+i,1), refh.offset(levelmin+i,2)); + LOGUSER(" size =(%5d,%5d,%5d)",refh.size(levelmin+i,0), + refh.size(levelmin+i,1), refh.size(levelmin+i,2)); + fine = new PaddedDensitySubGrid(refh.offset(levelmin+i,0), + refh.offset(levelmin+i,1), + refh.offset(levelmin+i,2), + refh.size(levelmin+i,0), + refh.size(levelmin+i,1), + refh.size(levelmin+i,2) ); + ///////////////////////////////////////////////////////////////////////// + + // load white noise for patch rand.load(*fine,levelmin+i); - convolution::perform( the_tf_kernel->fetch_kernel( levelmin+i, true ), reinterpret_cast( fine->get_data_ptr() ), shift ); + convolution::perform( the_tf_kernel->fetch_kernel( levelmin+i, true ), + reinterpret_cast( fine->get_data_ptr() ), shift ); if( i==1 ) fft_interpolate( *top, *fine, true ); else fft_interpolate( *coarse, *fine, false ); - - delta.add_patch( refh.offset(levelmin+i,0), refh.offset(levelmin+i,1), refh.offset(levelmin+i,2), - refh.size(levelmin+i,0), refh.size(levelmin+i,1), refh.size(levelmin+i,2) ); + delta.add_patch( refh.offset(levelmin+i,0), + refh.offset(levelmin+i,1), + refh.offset(levelmin+i,2), + refh.size(levelmin+i,0), + refh.size(levelmin+i,1), + refh.size(levelmin+i,2) ); fine->copy_unpad( *delta.get_grid(levelmin+i) ); - if( i==1 ) - delete top; - else - delete coarse; + if( i==1 ) delete top; + else delete coarse; coarse = fine; } @@ -466,7 +453,8 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type if( levelmax == levelmin ) { - std::cout << " - Performing noise convolution on level " << std::setw(2) << levelmax << " ..." << std::endl; + std::cout << " - Performing noise convolution on level " + << std::setw(2) << levelmax << " ..." << std::endl; LOGUSER("Performing noise convolution on level %3d...",levelmax); top = new DensityGrid( nbase, nbase, nbase ); @@ -474,7 +462,9 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type rand.load( *top, levelmin ); - convolution::perform( the_tf_kernel->fetch_kernel( levelmax ), reinterpret_cast( top->get_data_ptr() ), shift ); + convolution::perform( the_tf_kernel->fetch_kernel( levelmax ), + reinterpret_cast( top->get_data_ptr() ), + shift ); the_tf_kernel->deallocate(); delta.create_base_hierarchy(levelmin); @@ -496,8 +486,12 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type rand.load(*top,levelmin); } - fine = new PaddedDensitySubGrid( refh.offset(levelmin+i+1,0), refh.offset(levelmin+i+1,1), refh.offset(levelmin+i+1,2), - refh.size(levelmin+i+1,0), refh.size(levelmin+i+1,1), refh.size(levelmin+i+1,2) ); + fine = new PaddedDensitySubGrid( refh.offset(levelmin+i+1,0), + refh.offset(levelmin+i+1,1), + refh.offset(levelmin+i+1,2), + refh.size(levelmin+i+1,0), + refh.size(levelmin+i+1,1), + refh.size(levelmin+i+1,2) ); rand.load(*fine,levelmin+i+1); @@ -509,7 +503,8 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type /**********************************************************************************************************\ * multi-grid: top-level grid grids ..... \**********************************************************************************************************/ - std::cout << " - Performing noise convolution on level " << std::setw(2) << levelmin+i << " ..." << std::endl; + std::cout << " - Performing noise convolution on level " + << std::setw(2) << levelmin+i << " ..." << std::endl; LOGUSER("Performing noise convolution on level %3d",levelmin+i); LOGUSER("Creating base hierarchy..."); @@ -756,8 +751,8 @@ void coarsen_density( const refinement_hierarchy& rh, GridHierarchy& u ) } } - for( int i=rh.levelmax(); i>0; --i ) - mg_straight().restrict( *(u.get_grid(i)), *(u.get_grid(i-1)) ); + for( int i=rh.levelmax(); i>0; --i ) + mg_straight().restrict( *(u.get_grid(i)), *(u.get_grid(i-1)) ); } diff --git a/plugins/convex_hull.hh b/plugins/convex_hull.hh index 30a5227..eab7766 100644 --- a/plugins/convex_hull.hh +++ b/plugins/convex_hull.hh @@ -286,7 +286,7 @@ struct convex_hull{ for( j=i, l=0; l= 0 ) j=l; - int nt = omp_get_num_threads(); + int nt = omp_get_max_threads(); omp_set_num_threads( std::min(2,omp_get_max_threads()) ); #pragma omp parallel for diff --git a/plugins/output_gadget_tetmesh.cc b/plugins/output_gadget_tetmesh.cc index 8414bc5..feb4aa4 100644 --- a/plugins/output_gadget_tetmesh.cc +++ b/plugins/output_gadget_tetmesh.cc @@ -205,7 +205,7 @@ inline real_t get_cic( const grid_hierarchy & gh, int ilevel, real_t u, real_t v } -MyIDType compute_midpoint( const size_t *connect ) +MyIDType compute_midpoint( const MyIDType *connect ) { MyIDType lid1, lid2, newlid, lcoord1[3], lcoord2[3]; @@ -282,7 +282,7 @@ void split_lagrange_cube( size_t ip ) // insert mass carying "true" particles for( k=0; k #include #include +#include #include "output.hh" diff --git a/plugins/random_music.cc b/plugins/random_music.cc index 1745d33..3ec1689 100644 --- a/plugins/random_music.cc +++ b/plugins/random_music.cc @@ -5,10 +5,11 @@ public: explicit RNG_music( config_file& cf ) : RNG_plugin( cf ) { } + ~RNG_music() { } + bool is_multiscale() const - { - } + { return true; } }; diff --git a/plugins/region_convex_hull.cc b/plugins/region_convex_hull.cc index 6ae53dd..4a55c42 100644 --- a/plugins/region_convex_hull.cc +++ b/plugins/region_convex_hull.cc @@ -1,3 +1,13 @@ +/* + + region_convex_hull.cc - This file is part of MUSIC - + a code to generate multi-scale initial conditions + for cosmological simulations + + Copyright (C) 2010-13 Oliver Hahn + + */ + #include #include #include @@ -68,7 +78,7 @@ public: phull_->expand( sqrt(3.)*dx ); // output the center - float c[3] = { phull_->centroid_[0], phull_->centroid_[1], phull_->centroid_[2] }; + double c[3] = { phull_->centroid_[0], phull_->centroid_[1], phull_->centroid_[2] }; LOGINFO("Region center from convex hull centroid determined at\n\t (%f,%f,%f)",c[0],c[1],c[2]); //----------------------------------------------------------------- @@ -140,7 +150,7 @@ public: void get_center_unshifted( double *xc ) { double dx = 1.0/(1<centroid_[0], phull_->centroid_[1], phull_->centroid_[2] }; + double c[3] = { phull_->centroid_[0], phull_->centroid_[1], phull_->centroid_[2] }; xc[0] = c[0]+shift[0]*dx; xc[1] = c[1]+shift[1]*dx; xc[2] = c[2]+shift[2]*dx; diff --git a/plugins/region_ellipsoid.cc b/plugins/region_ellipsoid.cc index 33c4ff5..2838e96 100644 --- a/plugins/region_ellipsoid.cc +++ b/plugins/region_ellipsoid.cc @@ -1,4 +1,14 @@ #include +/* + + region_ellipsoid.cc - This file is part of MUSIC - + a code to generate multi-scale initial conditions + for cosmological simulations + + Copyright (C) 2010-13 Oliver Hahn + + */ + #include #include #include @@ -53,7 +63,7 @@ void Inverse_4x4( float *mat ) double det; /* determinant */ double dst[16]; -/* transpose matrix */ + /* transpose matrix */ for (int i = 0; i < 4; i++) { src[i] = mat[i*4]; @@ -397,9 +407,9 @@ public: { dist = (dist + 1.0) * detA13; - float q[3] = {x[0]-c[0],x[1]-c[1],x[2]-c[2]}; + T q[3] = {x[0]-c[0],x[1]-c[1],x[2]-c[2]}; - double r = 0.0; + T r = 0.0; for( int i=0; i<3; ++i ) for( int j=0; j<3; ++j ) r += q[i]*A[3*j+i]*q[j]; @@ -517,7 +527,7 @@ public: { std::vector pp; - for( int i=0; i<=levelmax_; ++i ) + for( unsigned i=0; i<=levelmax_; ++i ) pellip_.push_back( NULL ); @@ -537,12 +547,10 @@ public: padding_ = cf.getValue("setup","padding"); std::string point_file; - bool bfrom_file = true; - + if( cf.containsKey("setup", "region_point_file") ) { point_file = cf.getValue("setup","region_point_file"); - bfrom_file = true; point_reader pfr; pfr.read_points_from_file( point_file, vfac_, pp ); @@ -668,14 +676,14 @@ public: ~region_ellipsoid_plugin() { - for( int i=0; i<=levelmax_; ++i ) + for( unsigned i=0; i<=levelmax_; ++i ) if( pellip_[i] != NULL ) delete pellip_[i]; } void get_AABB( double *left, double *right, unsigned level ) { - if( level <= levelmin_ ) + if( level <= levelmin_ ) { left[0] = left[1] = left[2] = 0.0; right[0] = right[1] = right[2] = 1.0; diff --git a/region_generator.hh b/region_generator.hh index ebdc6e0..603578c 100644 --- a/region_generator.hh +++ b/region_generator.hh @@ -12,7 +12,7 @@ class region_generator_plugin{ public: config_file *pcf_; - int levelmin_, levelmax_; + unsigned levelmin_, levelmax_; public: region_generator_plugin( config_file& cf )