From 5d84806257cc343319712ef09eea673949235177 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Tue, 27 Aug 2013 21:50:51 -0700 Subject: [PATCH] fixes to the noise interpolation, improves stability of large-scale modes substantially --- Makefile | 2 +- main.cc | 2 +- random.cc | 90 ++++++++++++++++++++++++++++++++++++++++++------------- 3 files changed, 71 insertions(+), 23 deletions(-) diff --git a/Makefile b/Makefile index 36fcad5..275fad9 100644 --- a/Makefile +++ b/Makefile @@ -3,7 +3,7 @@ FFTW3 = yes MULTITHREADFFTW = yes SINGLEPRECISION = no -HAVEHDF5 = no +HAVEHDF5 = yes HAVEBOXLIB = no BOXLIB_HOME = ${HOME}/nyx_tot_sterben/BoxLib diff --git a/main.cc b/main.cc index 1436304..d1928a8 100644 --- a/main.cc +++ b/main.cc @@ -37,7 +37,7 @@ #include "transfer_function.hh" #define THE_CODE_NAME "music!" -#define THE_CODE_VERSION "1.3b" +#define THE_CODE_VERSION "1.4b" namespace music diff --git a/random.cc b/random.cc index 7a405c5..9e99171 100644 --- a/random.cc +++ b/random.cc @@ -499,7 +499,7 @@ random_numbers::random_numbers( random_numbers& rc, unsigned cubesize, lon { LOGINFO("Generating a constrained random number set with seed %ld\n using coarse mode replacement...",baseseed); - + assert(lx[0]%2==0 && lx[1]%2==0 && lx[2]%2==0); size_t nx=lx[0], ny=lx[1], nz=lx[2], nxc=lx[0]/2, nyc=lx[1]/2, nzc=lx[2]/2; @@ -575,25 +575,73 @@ random_numbers::random_numbers( random_numbers& rc, unsigned cubesize, lon #endif double fftnorm = 1.0/((double)nx*(double)ny*(double)nz); - - #pragma omp parallel for - for( int i=0; i<(int)nxc; i++ ) - for( int j=0; j<(int)nyc; j++ ) + double sqrt8 = sqrt(8.0); + // embedding + + // 0 0 +#pragma omp parallel for + for( int i=0; i<(int)nxc/2+1; i++ ) + for( int j=0; j<(int)nyc/2+1; j++ ) for( int k=0; k<(int)nzc/2+1; k++ ) - { - int ii(i),jj(j),kk(k); - - if( i > (int)nxc/2 ) ii += nx/2; - if( j > (int)nyc/2 ) jj += ny/2; - - size_t qc,qf; + { + int ii(i),jj(j),kk(k); + size_t qc,qf; qc = ((size_t)i*(size_t)nyc+(size_t)j)*(nzc/2+1)+(size_t)k; qf = ((size_t)ii*(size_t)ny+(size_t)jj)*(nz/2+1)+(size_t)kk; - - RE(cfine[qf]) = sqrt(8.0)*RE(ccoarse[qc]); - IM(cfine[qf]) = sqrt(8.0)*IM(ccoarse[qc]); - } - + + RE(cfine[qf]) = sqrt8*RE(ccoarse[qc]); + IM(cfine[qf]) = sqrt8*IM(ccoarse[qc]); + } + // 1 0 +#pragma omp parallel for + for( int i=nxc/2; i<(int)nxc; i++ ) + for( int j=0; j<(int)nyc/2+1; j++ ) + for( int k=0; k<(int)nzc/2+1; k++ ) + { + int ii(i+nx/2),jj(j),kk(k); + size_t qc,qf; + qc = ((size_t)i*(size_t)nyc+(size_t)j)*(nzc/2+1)+(size_t)k; + qf = ((size_t)ii*(size_t)ny+(size_t)jj)*(nz/2+1)+(size_t)kk; + + RE(cfine[qf]) = sqrt8*RE(ccoarse[qc]); + IM(cfine[qf]) = sqrt8*IM(ccoarse[qc]); + + if( i==(int)nxc/2 || j==(int)nyc/2 ) + IM(cfine[qf]) *= -1.0; + } + // 0 1 +#pragma omp parallel for + for( int i=0; i<(int)nxc/2+1; i++ ) + for( int j=nyc/2; j<(int)nyc; j++ ) + for( int k=0; k<(int)nzc/2+1; k++ ) + { + int ii(i),jj(j+ny/2),kk(k); + size_t qc,qf; + qc = ((size_t)i*(size_t)nyc+(size_t)j)*(nzc/2+1)+(size_t)k; + qf = ((size_t)ii*(size_t)ny+(size_t)jj)*(nz/2+1)+(size_t)kk; + + RE(cfine[qf]) = sqrt8*RE(ccoarse[qc]); + IM(cfine[qf]) = sqrt8*IM(ccoarse[qc]); + + if( i==(int)nxc/2 || j==(int)nyc/2 ) + IM(cfine[qf]) *= -1.0; + } + + // 1 1 +#pragma omp parallel for + for( int i=nxc/2; i<(int)nxc; i++ ) + for( int j=nyc/2; j<(int)nyc; j++ ) + for( int k=0; k<(int)nzc/2+1; k++ ) + { + int ii(i+nx/2),jj(j+ny/2),kk(k); + size_t qc,qf; + qc = ((size_t)i*(size_t)nyc+(size_t)j)*(nzc/2+1)+(size_t)k; + qf = ((size_t)ii*(size_t)ny+(size_t)jj)*(nz/2+1)+(size_t)kk; + + RE(cfine[qf]) = sqrt8*RE(ccoarse[qc]); + IM(cfine[qf]) = sqrt8*IM(ccoarse[qc]); + } + delete[] rcoarse; #pragma omp parallel for @@ -655,11 +703,11 @@ random_numbers::random_numbers( random_numbers& rc, unsigned cubesize, lon { LOGINFO("Generating a constrained random number set with seed %ld\n using Hoffman-Ribak constraints...", baseseed); - double fac = 1./sqrt(8.0); + double fac = 1.0/sqrt(8.0);//1./sqrt(8.0); - for( int i=x0[0],ii=x0[0]/2; ii