From 288749aea9cb4376aed02a22367a2d8078c7d0df Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Fri, 2 Jul 2010 16:17:14 -0700 Subject: [PATCH] minor fixes --- Makefile | 2 +- densities.cc | 15 +++------------ main.cc | 16 +++++++++++----- mg_interp.hh | 24 ++++++++++++------------ mg_solver.hh | 2 +- poisson.cc | 2 +- poisson.hh | 1 + 7 files changed, 30 insertions(+), 32 deletions(-) diff --git a/Makefile b/Makefile index 93a2de6..40afc5f 100644 --- a/Makefile +++ b/Makefile @@ -45,7 +45,7 @@ endif CFLAGS += $(OPT) TARGET = MUSIC OBJS = output.o transfer_function.o Numerics.o \ - convolution_kernel.o densities.o cosmology.o main.o \ + convolution_kernel.o densities.o cosmology.o poisson.o main.o \ $(patsubst plugins/%.cc,plugins/%.o,$(wildcard plugins/*.cc)) ############################################################################## diff --git a/densities.cc b/densities.cc index 88bc784..cb103c5 100644 --- a/densities.cc +++ b/densities.cc @@ -264,12 +264,8 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type #endif convolution::parameters conv_param; - //conv_param.boxlength = boxlength; conv_param.ptf = ptf; conv_param.pcf = &cf; - //conv_param.cosmo = cosmo; - - //... compute absolute grid offsets @@ -314,9 +310,6 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type //... if random numbers are not given for lower levels, obtain them by averaging if( lmingiven >= (int)levelmin ) { - //double x0[3] = { 0.0, 0.0, 0.0 }; - //double lx[3] = { 1.0, 1.0, 1.0 }; - randc[lmingiven] = new random_numbers( (unsigned)pow(2,lmingiven), 32, rngseeds[lmingiven], true );//, x0, lx ); for( int ilevel = lmingiven-1; ilevel >= (int)levelmin; --ilevel ){ @@ -332,8 +325,6 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type if( lmingiven < (int)levelmin ) { throw std::runtime_error("You provided a seed for a level below levelmin, this is not supported yet."); - //double x0[3] = { 0.0, 0.0, 0.0 }; - //double lx[3] = { 1.0, 1.0, 1.0 }; randc[lmingiven] = new random_numbers( (unsigned)pow(2,lmingiven), 32, rngseeds[lmingiven], true );//, x0, lx ); @@ -382,7 +373,7 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type random_numbers *rc; if( randc[levelmin] == NULL ) - rc = new random_numbers( nbase, 32, rngseeds[levelmin], true );//, x0, lx ); + rc = new random_numbers( nbase, 32, rngseeds[levelmin], true ); else rc = randc[levelmin]; @@ -760,7 +751,7 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type double sum; -#if 0 +#if 1 //... subtract the box mean.... this will otherwise add //... a constant curvature term to the potential sum = 0.0; @@ -779,7 +770,7 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type sum /= (nx*ny*nz); } - std::cout << " - Top grid mean density is off by " << sum << ", correcting..." << std::endl; + //std::cout << " - Top grid mean density is off by " << sum << ", correcting..." << std::endl; for( unsigned i=levelmin; i<=levelmax; ++i ) { diff --git a/main.cc b/main.cc index 244ecb7..dc0ff7d 100644 --- a/main.cc +++ b/main.cc @@ -48,7 +48,7 @@ #include "transfer_function.hh" #define THE_CODE_NAME "music!" -#define THE_CODE_VERSION "0.2a" +#define THE_CODE_VERSION "0.3a" namespace music @@ -86,15 +86,17 @@ void splash(void) void modify_grid_for_TF( const refinement_hierarchy& rh_full, refinement_hierarchy& rh_TF, config_file& cf ) { - unsigned lbase, lbaseTF, lmax; + unsigned lbase, lbaseTF, lmax, overlap; lbase = cf.getValue( "setup", "levelmin" ); lbaseTF = cf.getValueSafe( "setup", "levelminTF", lbase ); lmax = cf.getValue( "setup", "levelmax" ); + // TODO: add documentation about setup/overlap to manual + overlap = cf.getValueSafe( "setup", "overlap", 4 ); rh_TF = rh_full; - unsigned pad = 4; + unsigned pad = overlap; for( unsigned i=lbase+1; i<=lmax; ++i ) { @@ -287,6 +289,9 @@ int main (int argc, const char * argv[]) cosmo.nspect = cf.getValue( "cosmology", "nspec" ); cosmo.WDMg_x = cf.getValueSafe( "cosmology", "WDMg_x", 1.5 ); cosmo.WDMmass = cf.getValueSafe( "cosmology", "WDMmass", 0.0 ); + cosmo.dplus = 0.0; + cosmo.pnorm = 0.0; + cosmo.vfact = 0.0; //cosmo.Gamma = cf.getValueSafe( "cosmology", "Gamma", -1.0 ); @@ -567,10 +572,11 @@ int main (int argc, const char * argv[]) } the_output_plugin->finalize(); - }catch(...) + }catch(std::runtime_error& excp) { + std::cerr << " - " << excp.what() << std::endl; std::cerr << " - A fatal error occured. We need to exit...\n"; - bfatal; + bfatal = true; } //... clean up diff --git a/mg_interp.hh b/mg_interp.hh index 9afcad5..e81a349 100644 --- a/mg_interp.hh +++ b/mg_interp.hh @@ -692,7 +692,7 @@ struct interp_O5_fluxcorr ustar[q+2] = interp4( (*utop)(ixtop+p-1,iytop-2,iztop+q), (*utop)(ixtop+p-1,iytop-1,iztop+q), (*utop)(ixtop+p-1,iytop,iztop+q), (*utop)(ixtop+p-1,iytop+1,iztop+q), (*utop)(ixtop+p-1,iytop+2,iztop+q), fac*((double)j-0.5) ); - uhat[p] = interp4( ustar, fac*((double)k-0.5)-1.5 ); + uhat[p] = interp4( ustar, fac*((double)k-0.5) );//-1.5 ); } (*u)(ix,iy+j,iz+k) = interp4left( uhat[0], uhat[1], (*u)(ix+1,iy+j,iz+k), @@ -733,7 +733,7 @@ struct interp_O5_fluxcorr ustar[q+2] = interp4( (*utop)(ixtop+p,iytop-2,iztop+q), (*utop)(ixtop+p,iytop-1,iztop+q), (*utop)(ixtop+p,iytop,iztop+q), (*utop)(ixtop+p,iytop+1,iztop+q), (*utop)(ixtop+p,iytop+2,iztop+q), fac*((double)j-0.5) ); - uhat[p] = interp4( ustar, fac*((double)k-0.5)-1.5 ); + uhat[p] = interp4( ustar, fac*((double)k-0.5));//-1.5 ); } (*u)(ix,iy+j,iz+k) = interp4right( (*u)(ix-3,iy+j,iz+k), (*u)(ix-2,iy+j,iz+k), @@ -773,7 +773,7 @@ struct interp_O5_fluxcorr ustar[q+2] = interp4( (*utop)(ixtop-2,iytop+p-1,iztop+q), (*utop)(ixtop-1,iytop+p-1,iztop+q), (*utop)(ixtop,iytop+p-1,iztop+q), (*utop)(ixtop+1,iytop+p-1,iztop+q), (*utop)(ixtop+2,iytop+p-1,iztop+q), fac*((double)j-0.5) ); - uhat[p] = interp4( ustar, fac*((double)k-0.5)-1.5 ); + uhat[p] = interp4( ustar, fac*((double)k-0.5));//-1.5 ); } (*u)(ix+j,iy,iz+k) = interp4left( uhat[0], uhat[1], (*u)(ix+j,iy+1,iz+k), @@ -813,7 +813,7 @@ struct interp_O5_fluxcorr ustar[q+2] = interp4( (*utop)(ixtop-2,iytop+p,iztop+q), (*utop)(ixtop-1,iytop+p,iztop+q), (*utop)(ixtop,iytop+p,iztop+q), (*utop)(ixtop+1,iytop+p,iztop+q), (*utop)(ixtop+2,iytop+p,iztop+q), fac*((double)j-0.5) ); - uhat[p] = interp4( ustar, fac*((double)k-0.5)+1.5 ); + uhat[p] = interp4( ustar, fac*((double)k-0.5));//+1.5 ); } (*u)(ix+j,iy,iz+k) = interp4right( (*u)(ix+j,iy-3,iz+k), (*u)(ix+j,iy-2,iz+k), @@ -854,7 +854,7 @@ struct interp_O5_fluxcorr ustar[q+2] = interp4( (*utop)(ixtop-2,iytop+q,iztop+p-1), (*utop)(ixtop-1,iytop+q,iztop+p-1), (*utop)(ixtop,iytop+q,iztop+p-1), (*utop)(ixtop+1,iytop+q,iztop+p-1), (*utop)(ixtop+2,iytop+q,iztop+p-1), fac*((double)j-0.5) ); - uhat[p] = interp4( ustar, fac*((double)k-0.5)-1.5 ); + uhat[p] = interp4( ustar, fac*((double)k-0.5));//-1.5 ); } (*u)(ix+j,iy+k,iz) = interp4left( uhat[0], uhat[1], (*u)(ix+j,iy+k,iz+1), @@ -895,7 +895,7 @@ struct interp_O5_fluxcorr ustar[q+2] = interp4( (*utop)(ixtop-2,iytop+q,iztop+p), (*utop)(ixtop-1,iytop+q,iztop+p), (*utop)(ixtop,iytop+q,iztop+p), (*utop)(ixtop+1,iytop+q,iztop+p), (*utop)(ixtop+2,iytop+q,iztop+p), fac*((double)j-0.5) ); - uhat[p] = interp4( ustar, fac*((double)k-0.5)+1.5 ); + uhat[p] = interp4( ustar, fac*((double)k-0.5));//+1.5 ); } (*u)(ix+j,iy+k,iz) = interp4right( (*u)(ix+j,iy+k,iz-3), (*u)(ix+j,iy+k,iz-2), @@ -1000,7 +1000,7 @@ struct interp_O7_fluxcorr (*utop)(ixtop+p-2,iytop-1,iztop+q), (*utop)(ixtop+p-2,iytop,iztop+q), (*utop)(ixtop+p-2,iytop+1,iztop+q), (*utop)(ixtop+p-2,iytop+2,iztop+q), (*utop)(ixtop+p-2,iytop+3,iztop+q), fac*((double)j-0.5) ); - uhat[p] = interp6( ustar, fac*((double)k-0.5)-1.5 ); + uhat[p] = interp6( ustar, fac*((double)k-0.5));//-1.5 ); } (*u)(ix,iy+j,iz+k) = interp6left( uhat[0], uhat[1], uhat[2], (*u)(ix+1,iy+j,iz+k), @@ -1045,7 +1045,7 @@ struct interp_O7_fluxcorr (*utop)(ixtop+p,iytop-1,iztop+q), (*utop)(ixtop+p,iytop,iztop+q), (*utop)(ixtop+p,iytop+1,iztop+q), (*utop)(ixtop+p,iytop+2,iztop+q), (*utop)(ixtop+p,iytop+3,iztop+q), fac*((double)j-0.5) ); - uhat[p] = interp6( ustar, fac*((double)k-0.5)-1.5 ); + uhat[p] = interp6( ustar, fac*((double)k-0.5) );//-1.5 ); } (*u)(ix,iy+j,iz+k) = interp6right( (*u)(ix-4,iy+j,iz+k), (*u)(ix-3,iy+j,iz+k), (*u)(ix-2,iy+j,iz+k), @@ -1091,7 +1091,7 @@ struct interp_O7_fluxcorr (*utop)(ixtop-1,iytop+p-2,iztop+q), (*utop)(ixtop,iytop+p-2,iztop+q), (*utop)(ixtop+1,iytop+p-2,iztop+q), (*utop)(ixtop+2,iytop+p-2,iztop+q), (*utop)(ixtop+3,iytop+p-2,iztop+q), fac*((double)j-0.5) ); - uhat[p] = interp6( ustar, fac*((double)k-0.5)-1.5 ); + uhat[p] = interp6( ustar, fac*((double)k-0.5));//-1.5 ); } (*u)(ix+j,iy,iz+k) = interp6left( uhat[0], uhat[1], uhat[2], (*u)(ix+j,iy+1,iz+k), @@ -1136,7 +1136,7 @@ struct interp_O7_fluxcorr (*utop)(ixtop-1,iytop+p,iztop+q), (*utop)(ixtop,iytop+p,iztop+q), (*utop)(ixtop+1,iytop+p,iztop+q), (*utop)(ixtop+2,iytop+p,iztop+q), (*utop)(ixtop+3,iytop+p,iztop+q), fac*((double)j-0.5) ); - uhat[p] = interp6( ustar, fac*((double)k-0.5)+1.5 ); + uhat[p] = interp6( ustar, fac*((double)k-0.5));//+1.5 ); } (*u)(ix+j,iy,iz+k) = interp6right( (*u)(ix+j,iy-4,iz+k), (*u)(ix+j,iy-3,iz+k), (*u)(ix+j,iy-2,iz+k), @@ -1182,7 +1182,7 @@ struct interp_O7_fluxcorr (*utop)(ixtop-1,iytop+q,iztop+p-2), (*utop)(ixtop,iytop+q,iztop+p-2), (*utop)(ixtop+1,iytop+q,iztop+p-2), (*utop)(ixtop+2,iytop+q,iztop+p-2), (*utop)(ixtop+3,iytop+q,iztop+p-2), fac*((double)j-0.5) ); - uhat[p] = interp6( ustar, fac*((double)k-0.5)-1.5 ); + uhat[p] = interp6( ustar, fac*((double)k-0.5));//-1.5 ); } (*u)(ix+j,iy+k,iz) = interp6left( uhat[0], uhat[1], uhat[2], (*u)(ix+j,iy+k,iz+1), @@ -1227,7 +1227,7 @@ struct interp_O7_fluxcorr (*utop)(ixtop-1,iytop+q,iztop+p), (*utop)(ixtop,iytop+q,iztop+p), (*utop)(ixtop+1,iytop+q,iztop+p), (*utop)(ixtop+2,iytop+q,iztop+p), (*utop)(ixtop+3,iytop+q,iztop+p), fac*((double)j-0.5) ); - uhat[p] = interp6( ustar, fac*((double)k-0.5)+1.5 ); + uhat[p] = interp6( ustar, fac*((double)k-0.5));//+1.5 ); } (*u)(ix+j,iy+k,iz) = interp6right( (*u)(ix+j,iy+k,iz-4), (*u)(ix+j,iy+k,iz-3), (*u)(ix+j,iy+k,iz-2), diff --git a/mg_solver.hh b/mg_solver.hh index b21e09a..dcd44d8 100644 --- a/mg_solver.hh +++ b/mg_solver.hh @@ -114,7 +114,7 @@ template< class S, class I, class O, typename T > solver::solver( GridHierarchy& f, //const MeshvarBnd& ubnd, opt::smtype smoother, unsigned npresmooth, unsigned npostsmooth ) : m_scheme(), m_gridop(), m_npresmooth( npresmooth ), m_npostsmooth( npostsmooth ), -m_smoother( smoother ), m_ilevelmin( f.levelmin() ), m_pf( &f ), m_is_ini( true ) +m_smoother( smoother ), m_ilevelmin( f.levelmin() ), m_is_ini( true ), m_pf( &f ) { m_is_ini = true; diff --git a/poisson.cc b/poisson.cc index a5c6243..5d602fe 100644 --- a/poisson.cc +++ b/poisson.cc @@ -85,7 +85,7 @@ double multigrid_poisson_plugin::solve( grid_hierarchy& f, grid_hierarchy& u ) ps_presmooth = cf_.getValueSafe("poisson","pre_smooth",3); ps_postsmooth = cf_.getValueSafe("poisson","post_smooth",3); ps_smoother_name = cf_.getValueSafe("poisson","smoother","gs"); - order = cf_.getValueSafe( "setup", "laplace_order", 4 ); + order = cf_.getValueSafe( "poisson", "laplace_order", 4 ); multigrid::opt::smtype ps_smtype = multigrid::opt::sm_gauss_seidel; diff --git a/poisson.hh b/poisson.hh index 9c4d5cf..95e8ddb 100644 --- a/poisson.hh +++ b/poisson.hh @@ -113,6 +113,7 @@ protected: { double solve_O2( grid_hierarchy& f, grid_hierarchy& u ); double solve_O4( grid_hierarchy& f, grid_hierarchy& u ); + double solve_O6( grid_hierarchy& f, grid_hierarchy& u ); void gradient_O2( int dir, grid_hierarchy& u, grid_hierarchy& Du ); void gradient_O4( int dir, grid_hierarchy& u, grid_hierarchy& Du ); void gradient_O6( int dir, grid_hierarchy& u, grid_hierarchy& Du );