mirror/MUSIC
1
0
Fork 0
mirror of https://github.com/cosmo-sims/MUSIC.git synced 2024-09-19 17:03:46 +02:00
Code Releases Activity

Merged in panphasia (pull request #3)

Browse source 
Panphasia
This commit is contained in:
Oliver Hahn 2023-02-14 00:12:25 +00:00
commit 5304b8141a
19 changed files with 7268 additions and 2831 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

20
CMakeLists.txt
View file

@ -71,6 +71,22 @@ file( GLOB PLUGINS
${PROJECT_SOURCE_DIR}/src/plugins/*.cc
)
# PANPHASIA
option(ENABLE_PANPHASIA "Enable PANPHASIA random number generator" ON)
if(ENABLE_PANPHASIA)
enable_language(Fortran)
if ("${CMAKE_Fortran_COMPILER_ID}" MATCHES "Intel")
set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -132 -implicit-none")
elseif("${CMAKE_Fortran_COMPILER_ID}" MATCHES "GNU")
set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -ffixed-line-length-132 -fimplicit-none")
endif()
list (APPEND SOURCES
${PROJECT_SOURCE_DIR}/ext/panphasia/panphasia_routines.f
${PROJECT_SOURCE_DIR}/ext/panphasia/generic_lecuyer.f90
)
# target_include_directories(${PRGNAME} PRIVATE ${PROJECT_SOURCE_DIR}/external/panphasia_ho)
endif(ENABLE_PANPHASIA)
add_executable(${PRGNAME} ${SOURCES} ${PLUGINS})
set_target_properties(${PRGNAME} PROPERTIES CXX_STANDARD 11)
@ -116,6 +132,10 @@ if(TIRPC_FOUND)
target_compile_options(${PRGNAME} PRIVATE "-DHAVE_TIRPC")
endif(TIRPC_FOUND)
if(ENABLE_PANPHASIA)
target_compile_options(${PRGNAME} PRIVATE "-DHAVE_PANPHASIA")
endif(ENABLE_PANPHASIA)
target_link_libraries(${PRGNAME} ${FFTW3_LIBRARIES})
target_include_directories(${PRGNAME} PRIVATE ${FFTW3_INCLUDE_DIRS})

683
ext/panphasia/generic_lecuyer.f90 Normal file
View file

@ -0,0 +1,683 @@
!=====================================================================================c
!
! The code below was written by: Stephen Booth
! Edinburgh Parallel Computing Centre
! The University of Edinburgh
! JCMB
! Mayfield Road
! Edinburgh EH9 3JZ
! United Kingdom
!
! This file is part of the software made public in
! Jenkins and Booth 2013 - arXiv:1306.XXXX
!
! The software computes the Panphasia Gaussian white noise field
! realisation described in detail in Jenkins 2013 - arXiv:1306.XXXX
!
!
!
! This software is free, subject to a agreeing licence conditions:
!
!
! (i) you will publish the phase descriptors and reference Jenkins (13)
! for any new simulations that use Panphasia phases. You will pass on this
! condition to others for any software or data you make available publically
! or privately that makes use of Panphasia.
!
! (ii) that you will ensure any publications using results derived from Panphasia
! will be submitted as a final version to arXiv prior to or coincident with
! publication in a journal.
!
!
! (iii) that you report any bugs in this software as soon as confirmed to
! A.R.Jenkins@durham.ac.uk
!
! (iv) that you understand that this software comes with no warranty and that is
! your responsibility to ensure that it is suitable for the purpose that
! you intend.
!
!=====================================================================================c
!{{{Rand_base (define kind types)
MODULE Rand_base
! This module just declares the base types
! we may have to edit this to match to the target machine
! we really need a power of 2 selected int kind in fortran-95 we could
! do this with a PURE function I think.
!
! 10 decimal digits will hold 2^31
!
INTEGER, PARAMETER :: Sint = SELECTED_INT_KIND(9)
! INTEGER, PARAMETER :: Sint = SELECTED_INT_KIND(10)
! INTEGER, PARAMETER :: Sint = 4
!
! 18-19 decimal digits will hold 2^63
! but all 19 digit numbers require 2^65 :-(
!
INTEGER, PARAMETER :: Dint = SELECTED_INT_KIND(17)
! INTEGER, PARAMETER :: Dint = SELECTED_INT_KIND(18)
! INTEGER, PARAMETER :: Dint = 8
! type for index counters must hold Nstore
INTEGER, PARAMETER :: Ctype = SELECTED_INT_KIND(3)
END MODULE Rand_base
!}}}
!{{{Rand_int (random integers mod 2^31-1)
MODULE Rand_int
USE Rand_base
IMPLICIT NONE
! The general approach of this module is two have
! two types Sint and Dint
!
! Sint should have at least 31 bits
! dint shouldhave at least 63
!{{{constants
INTEGER(KIND=Ctype), PARAMETER :: Nstate=5_Ctype
INTEGER(KIND=Ctype), PRIVATE, PARAMETER :: Nbatch=128_Ctype
INTEGER(KIND=Ctype), PRIVATE, PARAMETER :: Nstore=Nstate+Nbatch
INTEGER(KIND=Sint), PRIVATE, PARAMETER :: M = 2147483647_Sint
INTEGER(KIND=Dint), PRIVATE, PARAMETER :: Mask = 2147483647_Dint
INTEGER(KIND=Dint), PRIVATE, PARAMETER :: A1 = 107374182_Dint
INTEGER(KIND=Dint), PRIVATE, PARAMETER :: A5 = 104480_Dint
LOGICAL, PARAMETER :: Can_step_int=.TRUE.
LOGICAL, PARAMETER :: Can_reverse_int=.TRUE.
!}}}
!{{{Types
!
! This type holds the state of the generator
!
!{{{TYPE RAND_state
TYPE RAND_state
PRIVATE
INTEGER(KIND=Sint) :: state(Nstore)
! do we need to re-fill state table this is reset when we initialise state.
LOGICAL :: need_fill
! position of the next state variable to output
INTEGER(KIND=Ctype) :: pos
END TYPE RAND_state
!}}}
!
! This type defines the offset type used for stepping.
!
!{{{TYPE RAND_offset
TYPE RAND_offset
PRIVATE
INTEGER(KIND=Sint) :: poly(Nstate)
END TYPE RAND_offset
!}}}
!}}}
!{{{interface and overloads
!
! Allow automatic conversion between integers and offsets
!
INTERFACE ASSIGNMENT(=)
MODULE PROCEDURE Rand_set_offset
MODULE PROCEDURE Rand_load
MODULE PROCEDURE Rand_save
MODULE PROCEDURE Rand_seed
END INTERFACE
INTERFACE OPERATOR(+)
MODULE PROCEDURE Rand_add_offset
END INTERFACE
INTERFACE OPERATOR(*)
MODULE PROCEDURE Rand_mul_offset
END INTERFACE
!
! overload + as the boost/stepping operator
!
INTERFACE OPERATOR(+)
MODULE PROCEDURE Rand_step
MODULE PROCEDURE Rand_boost
END INTERFACE
!}}}
!{{{PUBLIC/PRIVATE
PRIVATE reduce,mod_saxpy,mod_sdot,p_saxpy,p_sdot,poly_mult
PRIVATE poly_square, poly_power
PRIVATE fill_state, repack_state
PUBLIC Rand_sint, Rand_sint_vec
PUBLIC Rand_save, Rand_load
PUBLIC Rand_set_offset, Rand_add_offset, Rand_mul_offset
PUBLIC Rand_step, Rand_boost, Rand_seed
!}}}
CONTAINS
!{{{Internals
!{{{RECURSIVE FUNCTION reduce(A)
RECURSIVE FUNCTION reduce(A)
!
! Take A Dint and reduce to Sint MOD M
!
INTEGER(KIND=Dint), INTENT(IN) :: A
INTEGER(KIND=Sint) reduce
INTEGER(KIND=Dint) tmp
tmp = A
DO WHILE( ISHFT(tmp, -31) .GT. 0 )
tmp = IAND(tmp,Mask) + ISHFT(tmp, -31)
END DO
IF( tmp .GE. M ) THEN
reduce = tmp - M
ELSE
reduce = tmp
END IF
END FUNCTION reduce
!}}}
!{{{RECURSIVE SUBROUTINE fill_state(x)
RECURSIVE SUBROUTINE fill_state(x)
TYPE(RAND_state), INTENT(INOUT) :: x
INTEGER(KIND=Ctype) i
INTRINSIC IAND, ISHFT
INTEGER(KIND=Dint) tmp
DO i=Nstate+1,Nstore
tmp = (x%state(i-5) * A5) + (x%state(i-1)*A1)
!
! now reduce down to mod M efficiently
! really hope the compiler in-lines this
!
! x%state(i) = reduce(tmp)
DO WHILE( ISHFT(tmp, -31) .GT. 0 )
tmp = IAND(tmp,Mask) + ISHFT(tmp, -31)
END DO
IF( tmp .GE. M ) THEN
x%state(i) = tmp - M
ELSE
x%state(i) = tmp
END IF
END DO
x%need_fill = .FALSE.
END SUBROUTINE fill_state
!}}}
!{{{RECURSIVE SUBROUTINE repack_state(x)
RECURSIVE SUBROUTINE repack_state(x)
TYPE(RAND_state), INTENT(INOUT) :: x
INTEGER(KIND=Ctype) i
DO i=1,Nstate
x%state(i) = x%state(i+x%pos-(Nstate+1))
END DO
x%pos = Nstate + 1
x%need_fill = .TRUE.
END SUBROUTINE repack_state
!}}}
!{{{RECURSIVE SUBROUTINE mod_saxpy(y,a,x)
RECURSIVE SUBROUTINE mod_saxpy(y,a,x)
INTEGER(KIND=Ctype) i
INTEGER(KIND=Sint) y(Nstate)
INTEGER(KIND=Sint) a
INTEGER(KIND=Sint) x(Nstate)
INTEGER(KIND=Dint) tx,ty,ta
IF( a .EQ. 0_Sint ) RETURN
! We use KIND=Dint temporaries here to ensure
! that we don't overflow in the expression
ta = a
DO i=1,Nstate
ty=y(i)
tx=x(i)
y(i) = reduce(ty + ta * tx)
END DO
END SUBROUTINE
!}}}
!{{{RECURSIVE SUBROUTINE mod_sdot(res,x,y)
RECURSIVE SUBROUTINE mod_sdot(res,x,y)
INTEGER(KIND=Sint), INTENT(OUT) :: res
INTEGER(KIND=Sint), INTENT(IN) :: x(Nstate) , y(Nstate)
INTEGER(KIND=Dint) dx, dy, dtmp
INTEGER(KIND=Sint) tmp
INTEGER(KIND=Ctype) i
tmp = 0
DO i=1,Nstate
dx = x(i)
dy = y(i)
dtmp = tmp
tmp = reduce(dtmp + dx * dy)
END DO
res = tmp
END SUBROUTINE
!}}}
!{{{RECURSIVE SUBROUTINE p_saxpy(y,a)
RECURSIVE SUBROUTINE p_saxpy(y,a)
! Calculates mod_saxpy(y,a,P)
INTEGER(KIND=Sint), INTENT(INOUT) :: y(Nstate)
INTEGER(KIND=Sint), INTENT(IN) :: a
INTEGER(KIND=Dint) tmp, dy, da
dy = y(1)
da = a
tmp = dy + da*A5
y(1) = reduce(tmp)
dy = y(5)
da = a
tmp = dy + da*A1
y(5) = reduce(tmp)
END SUBROUTINE
!}}}
!{{{RECURSIVE SUBROUTINE p_sdot(res,n,x)
RECURSIVE SUBROUTINE p_sdot(res,x)
INTEGER(KIND=Sint), INTENT(OUT) :: res
INTEGER(KIND=Sint), INTENT(IN) :: x(Nstate)
INTEGER(KIND=Dint) dx1, dx5, dtmp
dx1 = x(1)
dx5 = x(5)
dtmp = A1*dx5 + A5*dx1
res = reduce(dtmp)
END SUBROUTINE
!}}}
!{{{RECURSIVE SUBROUTINE poly_mult(a,b)
RECURSIVE SUBROUTINE poly_mult(a,b)
INTEGER(KIND=Sint), INTENT(INOUT) :: a(Nstate)
INTEGER(KIND=Sint), INTENT(IN) :: b(Nstate)
INTEGER(KIND=Sint) tmp((2*Nstate) - 1)
INTEGER(KIND=Ctype) i
tmp = 0_Sint
DO i=1,Nstate
CALL mod_saxpy(tmp(i:Nstate+i-1),a(i), b)
END DO
DO i=(2*Nstate)-1, Nstate+1, -1
CALL P_SAXPY(tmp(i-Nstate:i-1),tmp(i))
END DO
a = tmp(1:Nstate)
END SUBROUTINE
!}}}
!{{{RECURSIVE SUBROUTINE poly_square(a)
RECURSIVE SUBROUTINE poly_square(a)
INTEGER(KIND=Sint), INTENT(INOUT) :: a(Nstate)
INTEGER(KIND=Sint) tmp((2*Nstate) - 1)
INTEGER(KIND=Ctype) i
tmp = 0_Sint
DO i=1,Nstate
CALL mod_saxpy(tmp(i:Nstate+i-1),a(i), a)
END DO
DO i=(2*Nstate)-1, Nstate+1, -1
CALL P_SAXPY(tmp(i-Nstate:i-1),tmp(i))
END DO
a = tmp(1:Nstate)
END SUBROUTINE
!}}}
!{{{RECURSIVE SUBROUTINE poly_power(poly,n)
RECURSIVE SUBROUTINE poly_power(poly,n)
INTEGER(KIND=Sint), INTENT(INOUT) :: poly(Nstate)
INTEGER, INTENT(IN) :: n
INTEGER nn
INTEGER(KIND=Sint) x(Nstate), out(Nstate)
IF( n .EQ. 0 )THEN
poly = 0_Sint
poly(1) = 1_Sint
RETURN
ELSE IF( n .LT. 0 )THEN
poly = 0_Sint
RETURN
END IF
out = 0_sint
out(1) = 1_Sint
x = poly
nn = n
DO WHILE( nn .GT. 0 )
IF( MOD(nn,2) .EQ. 1 )THEN
call poly_mult(out,x)
END IF
nn = nn/2
IF( nn .GT. 0 )THEN
call poly_square(x)
END IF
END DO
poly = out
END SUBROUTINE poly_power
!}}}
!}}}
!{{{RECURSIVE SUBROUTINE Rand_seed( state, n )
RECURSIVE SUBROUTINE Rand_seed( state, n )
TYPE(Rand_state), INTENT(OUT) :: state
INTEGER, INTENT(IN) :: n
! initialise the genrator using a single integer
! fist initialise to an arbitrary state then boost by a multiple
! of a long distance
!
! state is moved forward by P^n steps
! we want this to be ok for seperating parallel sequences on MPP machines
! P is taken as a prime number as this should prevent strong correlations
! when the generators are operated in tight lockstep.
! equivalent points on different processors will also be related by a
! primative polynomial
! P is 2^48-59
TYPE(Rand_state) tmp
TYPE(Rand_offset), PARAMETER :: P = &
Rand_offset( (/ 1509238949_Sint ,2146167999_Sint ,1539340803_Sint , &
1041407428_Sint ,666274987_Sint /) )
CALL Rand_load( tmp, (/ 5, 4, 3, 2, 1 /) )
state = Rand_boost( tmp, Rand_mul_offset(P, n ))
END SUBROUTINE Rand_seed
!}}}
!{{{RECURSIVE SUBROUTINE Rand_load( state, input )
RECURSIVE SUBROUTINE Rand_load( state, input )
TYPE(RAND_state), INTENT(OUT) :: state
INTEGER, INTENT(IN) :: input(Nstate)
INTEGER(KIND=Ctype) i
state%state = 0_Sint
DO i=1,Nstate
state%state(i) = MOD(INT(input(i),KIND=Sint),M)
END DO
state%need_fill = .TRUE.
state%pos = Nstate + 1
END SUBROUTINE Rand_load
!}}}
!{{{RECURSIVE SUBROUTINE Rand_save( save_vec,state )
RECURSIVE SUBROUTINE Rand_save( save_vec, x )
INTEGER, INTENT(OUT) :: save_vec(Nstate)
TYPE(RAND_state), INTENT(IN) :: x
INTEGER(KIND=Ctype) i
DO i=1,Nstate
save_vec(i) = x%state(x%pos-(Nstate+1) + i)
END DO
END SUBROUTINE Rand_save
!}}}
!{{{RECURSIVE SUBROUTINE Rand_set_offset( offset, n )
RECURSIVE SUBROUTINE Rand_set_offset( offset, n )
TYPE(Rand_offset), INTENT(OUT) :: offset
INTEGER, INTENT(IN) :: n
offset%poly = 0_Sint
IF ( n .GE. 0 ) THEN
offset%poly(2) = 1_Sint
call poly_power(offset%poly,n)
ELSE
!
! This is X^-1
!
offset%poly(4) = 858869107_Sint
offset%poly(5) = 1840344978_Sint
call poly_power(offset%poly,-n)
END IF
END SUBROUTINE Rand_set_offset
!}}}
!{{{TYPE(Rand_offset) RECURSIVE FUNCTION Rand_add_offset( a, b )
TYPE(Rand_offset) RECURSIVE FUNCTION Rand_add_offset( a, b )
TYPE(Rand_offset), INTENT(IN) :: a, b
Rand_add_offset = a
CALL poly_mult(Rand_add_offset%poly,b%poly)
RETURN
END FUNCTION Rand_add_offset
!}}}
!{{{TYPE(Rand_offset) RECURSIVE FUNCTION Rand_mul_offset( a, n )
TYPE(Rand_offset) RECURSIVE FUNCTION Rand_mul_offset( a, n )
TYPE(Rand_offset), INTENT(IN) :: a
INTEGER, INTENT(IN) :: n
Rand_mul_offset = a
CALL poly_power(Rand_mul_offset%poly,n)
RETURN
END FUNCTION Rand_mul_offset
!}}}
!{{{RECURSIVE FUNCTION Rand_boost(x, offset)
RECURSIVE FUNCTION Rand_boost(x, offset)
TYPE(Rand_state) Rand_boost
TYPE(Rand_state), INTENT(IN) :: x
TYPE(Rand_offset), INTENT(IN) :: offset
INTEGER(KIND=Sint) tmp(2*Nstate-1), res(Nstate)
INTEGER(KIND=Ctype) i
DO i=1,Nstate
tmp(i) = x%state(x%pos-(Nstate+1) + i)
END DO
tmp(Nstate+1:) = 0_Sint
DO i=1,Nstate-1
call P_SDOT(tmp(i+Nstate),tmp(i:Nstate+i-1))
END DO
DO i=1,Nstate
call mod_sdot(res(i),offset%poly,tmp(i:Nstate+i-1))
END DO
Rand_boost%state = 0_Sint
DO i=1,Nstate
Rand_boost%state(i) = res(i)
END DO
Rand_boost%need_fill = .TRUE.
Rand_boost%pos = Nstate + 1
END FUNCTION Rand_boost
!}}}
!{{{RECURSIVE FUNCTION Rand_step(x, n)
RECURSIVE FUNCTION Rand_step(x, n)
TYPE(Rand_state) Rand_step
TYPE(RAND_state), INTENT(IN) :: x
INTEGER, INTENT(IN) :: n
TYPE(Rand_offset) tmp
CALL Rand_set_offset(tmp,n)
Rand_step=Rand_boost(x,tmp)
END FUNCTION
!}}}
!{{{RECURSIVE FUNCTION Rand_sint(x)
RECURSIVE FUNCTION Rand_sint(x)
TYPE(RAND_state), INTENT(INOUT) :: x
INTEGER(KIND=Sint) Rand_sint
IF( x%pos .GT. Nstore )THEN
CALL repack_state(x)
END IF
IF( x%need_fill ) CALL fill_state(x)
Rand_sint = x%state(x%pos)
x%pos = x%pos + 1
RETURN
END FUNCTION Rand_sint
!}}}
!{{{RECURSIVE SUBROUTINE Rand_sint_vec(iv,x)
RECURSIVE SUBROUTINE Rand_sint_vec(iv,x)
INTEGER(KIND=Sint), INTENT(OUT) :: iv(:)
TYPE(RAND_state), INTENT(INOUT) :: x
INTEGER left,start, chunk, i
start=1
left=SIZE(iv)
DO WHILE( left .GT. 0 )
IF( x%pos .GT. Nstore )THEN
CALL repack_state(x)
END IF
IF( x%need_fill ) CALL fill_state(x)
chunk = MIN(left,Nstore-x%pos+1)
DO i=0,chunk-1
iv(start+i) = x%state(x%pos+i)
END DO
start = start + chunk
x%pos = x%pos + chunk
left = left - chunk
END DO
RETURN
END SUBROUTINE Rand_sint_vec
!}}}
END MODULE Rand_int
!}}}
!{{{Rand (use Rand_int to make random reals)
MODULE Rand
USE Rand_int
IMPLICIT NONE
!{{{Parameters
INTEGER, PARAMETER :: RAND_kind1 = SELECTED_REAL_KIND(10)
INTEGER, PARAMETER :: RAND_kind2 = SELECTED_REAL_KIND(6)
INTEGER, PARAMETER, PRIVATE :: Max_block=100
INTEGER(KIND=Sint), PRIVATE, PARAMETER :: M = 2147483647
REAL(KIND=RAND_kind1), PRIVATE, PARAMETER :: INVMP1_1 = ( 1.0_RAND_kind1 / 2147483647.0_RAND_kind1 )
REAL(KIND=RAND_kind2), PRIVATE, PARAMETER :: INVMP1_2 = ( 1.0_RAND_kind2 / 2147483647.0_RAND_kind2 )
LOGICAL, PARAMETER :: Can_step = Can_step_int
LOGICAL, PARAMETER :: Can_reverse = Can_reverse_int
!}}}
PUBLIC Rand_real
INTERFACE Rand_real
MODULE PROCEDURE Rand_real1
MODULE PROCEDURE Rand_real2
MODULE PROCEDURE Rand_real_vec1
MODULE PROCEDURE Rand_real_vec2
END INTERFACE
CONTAINS
!{{{RECURSIVE SUBROUTINE Rand_real1(y,x)
RECURSIVE SUBROUTINE Rand_real1(y,x)
REAL(KIND=RAND_kind1), INTENT(OUT) :: y
TYPE(RAND_state), INTENT(INOUT) :: x
INTEGER(KIND=Sint) Z
Z = Rand_sint(x)
IF (Z .EQ. 0) Z = M
y = ((Z-0.5d0)*INVMP1_1)
RETURN
END SUBROUTINE Rand_real1
!}}}
!{{{RECURSIVE SUBROUTINE Rand_real2(y,x)
RECURSIVE SUBROUTINE Rand_real2(y,x)
REAL(KIND=RAND_kind2), INTENT(OUT) :: y
TYPE(RAND_state), INTENT(INOUT) :: x
INTEGER(KIND=Sint) Z
Z = Rand_sint(x)
IF (Z .EQ. 0) Z = M
y = ((Z-0.5d0)*INVMP1_1) ! generate in double and truncate.
RETURN
END SUBROUTINE Rand_real2
!}}}
!{{{RECURSIVE SUBROUTINE Rand_real_vec1(rv,x)
RECURSIVE SUBROUTINE Rand_real_vec1(rv,x)
TYPE(RAND_state), INTENT(INOUT) :: x
REAL(KIND=RAND_kind1) rv(:)
INTEGER left,start, chunk, i
INTEGER(KIND=Sint) Z
INTEGER(KIND=Sint) temp(MIN(SIZE(rv),Max_block))
start=0
left=SIZE(rv)
DO WHILE( left .GT. 0 )
chunk = MIN(left,Max_block)
CALL Rand_sint_vec(temp(1:chunk),x)
DO i=1,chunk
Z = temp(i)
IF (Z .EQ. 0) Z = M
rv(start+i) = (Z-0.5d0)*INVMP1_1
END DO
start = start + chunk
left = left - chunk
END DO
RETURN
END SUBROUTINE Rand_real_vec1
!}}}
!{{{RECURSIVE SUBROUTINE Rand_real_vec2(rv,x)
RECURSIVE SUBROUTINE Rand_real_vec2(rv,x)
TYPE(RAND_state), INTENT(INOUT) :: x
REAL(KIND=RAND_kind2) rv(:)
INTEGER left,start, chunk, i
INTEGER(KIND=Sint) Z
INTEGER(KIND=Sint) temp(MIN(SIZE(rv),Max_block))
start=0
left=SIZE(rv)
DO WHILE( left .GT. 0 )
chunk = MIN(left,Max_block)
CALL Rand_sint_vec(temp(1:chunk),x)
DO i=1,chunk
Z = temp(i)
IF (Z .EQ. 0) Z = M
rv(start+i) = (Z-0.5d0)*INVMP1_2
END DO
start = start + chunk
left = left - chunk
END DO
RETURN
END SUBROUTINE Rand_real_vec2
!}}}
END MODULE Rand
!}}}
!{{{test program
! PROGRAM test_random
! use Rand
! TYPE(RAND_state) x
! REAL y
! CALL Rand_load(x,(/5,4,3,2,1/))
! DO I=0,10
! CALL Rand_real(y,x)
! WRITE(*,10) I,y
! END DO
!
!10 FORMAT(I10,E25.16)
!
! END
! 0 0.5024326127022505E-01
! 1 0.8260946767404675E-01
! 2 0.2123264316469431E-01
! 3 0.6926658791489899E+00
! 4 0.2076155943796039E+00
! 5 0.4327449947595596E-01
! 6 0.2204052871093154E-01
! 7 0.1288446951657534E+00
! 8 0.4859915426932275E+00
! 9 0.5721384193748236E-01
! 10 0.7996825082227588E+00
!
!}}}

3334
ext/panphasia/panphasia_routines.f Normal file
View file

File diff suppressed because it is too large Load diff

3
ics_example.conf
View file

@ -12,8 +12,7 @@ align_top = no
baryons = no
use_2LPT = no
use_LLA = no
periodic_TF = yes
zero_zoom_velocity = no
[cosmology]
Omega_m = 0.276

58
src/convolution_kernel.cc
View file

@ -1,10 +1,9 @@
/*
convolution_kernel.cc - This file is part of MUSIC -
a code to generate multi-scale initial conditions
for cosmological simulations
a code to generate multi-scale initial conditions for cosmological simulations
Copyright (C) 2010-19 Oliver Hahn
Copyright (C) 2010-23 Oliver Hahn
*/
@ -35,7 +34,8 @@ void perform(kernel *pk, void *pd, bool shift, bool fix, bool flip)
double fftnormp = 1.0/sqrt((double)cparam_.nx * (double)cparam_.ny * (double)cparam_.nz);
double fftnorm = pow(2.0 * M_PI, 1.5) / sqrt(cparam_.lx * cparam_.ly * cparam_.lz) * fftnormp;
fftw_complex *cdata, *ckernel;
fftw_complex *cdata;
[[maybe_unused]] fftw_complex *ckernel;
fftw_real *data;
data = reinterpret_cast<fftw_real *>(pd);
@ -96,51 +96,8 @@ void perform(kernel *pk, void *pd, bool shift, bool fix, bool flip)
std::complex<double> dcmode(RE(cdata[0]), IM(cdata[0]));
if (!pk->is_ksampled())
{
#pragma omp parallel for
for (int i = 0; i < cparam_.nx; ++i)
for (int j = 0; j < cparam_.ny; ++j)
for (int k = 0; k < cparam_.nz / 2 + 1; ++k)
{
size_t ii = (size_t)(i * cparam_.ny + j) * (size_t)(cparam_.nz / 2 + 1) + (size_t)k;
double kx, ky, kz;
kx = (double)i;
ky = (double)j;
kz = (double)k;
if (kx > cparam_.nx / 2)
kx -= cparam_.nx;
if (ky > cparam_.ny / 2)
ky -= cparam_.ny;
double arg = (kx + ky + kz) * dstag;
std::complex<double> carg(cos(arg), sin(arg));
std::complex<double>
ccdata(RE(cdata[ii]), IM(cdata[ii])),
cckernel(RE(ckernel[ii]), IM(ckernel[ii]));
if( fix ){
ccdata = ccdata / std::abs(ccdata);
}
if( flip ){
ccdata = -ccdata;
}
ccdata = ccdata * cckernel * fftnorm * carg;
RE(cdata[ii]) = ccdata.real();
IM(cdata[ii]) = ccdata.imag();
}
}
else
{
#pragma omp parallel
#pragma omp parallel
{
const size_t veclen = cparam_.nz / 2 + 1;
@ -149,7 +106,7 @@ void perform(kernel *pk, void *pd, bool shift, bool fix, bool flip)
double *Tkvec = new double[veclen];
double *argvec = new double[veclen];
#pragma omp for
#pragma omp for
for (int i = 0; i < cparam_.nx; ++i)
for (int j = 0; j < cparam_.ny; ++j)
{
@ -202,7 +159,7 @@ void perform(kernel *pk, void *pd, bool shift, bool fix, bool flip)
// we now set the correct DC mode below...
RE(cdata[0]) = 0.0;
IM(cdata[0]) = 0.0;
}
LOGUSER("Performing backward FFT...");
@ -229,7 +186,6 @@ void perform(kernel *pk, void *pd, bool shift, bool fix, bool flip)
#endif
// set the DC mode here to avoid a possible truncation error in single precision
if (pk->is_ksampled())
{
size_t nelem = (size_t)cparam_.nx * (size_t)cparam_.ny * (size_t)cparam_.nz;
real_t mean = dcmode.real() * fftnorm / (real_t)nelem;

3
src/convolution_kernel.hh
View file

@ -4,7 +4,7 @@
a code to generate multi-scale initial conditions
for cosmological simulations
Copyright (C) 2010-19 Oliver Hahn
Copyright (C) 2010-23 Oliver Hahn
*/
@ -16,6 +16,7 @@
#include "config_file.hh"
#include "densities.hh"
#include "density_grid.hh"
#include "transfer_function.hh"
#define ACC_RF(i, j, k) (((((size_t)(i) + nx) % nx) * ny + (((size_t)(j) + ny) % ny)) * 2 * (nz / 2 + 1) + (((size_t)(k) + nz) % nz))

45
src/densities.cc
View file

@ -11,6 +11,7 @@
#include <cstring>
#include "densities.hh"
#include "random.hh"
#include "convolution_kernel.hh"
//TODO: this should be a larger number by default, just to maintain consistency with old default
@ -113,8 +114,13 @@ void fft_coarsen(m1 &v, m2 &V)
val_fine *= val_phas * fftnorm / 8.0;
if( i!=(int)nxF/2 && j!=(int)nyF/2 && k!=(int)nzF/2 ){
RE(ccoarse[qc]) = val_fine.real();
IM(ccoarse[qc]) = val_fine.imag();
}else{
RE(ccoarse[qc]) = 0.0;//val_fine.real();
IM(ccoarse[qc]) = 0.0;//val_fine.imag();
}
}
delete[] rfine;
@ -335,7 +341,7 @@ void fft_interpolate(m1 &V, m2 &v, bool from_basegrid = false)
/*******************************************************************************************/
void GenerateDensityUnigrid(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, noise_generator &rand, grid_hierarchy &delta, bool smooth, bool shift)
{
unsigned levelmin, levelmax, levelminPoisson;
@ -416,7 +422,7 @@ void GenerateDensityUnigrid(config_file &cf, transfer_function *ptf, tf_type typ
/*******************************************************************************************/
void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type type,
refinement_hierarchy &refh, rand_gen &rand,
refinement_hierarchy &refh, noise_generator &rand,
grid_hierarchy &delta, bool smooth, bool shift)
{
unsigned levelmin, levelmax, levelminPoisson;
@ -439,6 +445,7 @@ void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type t
bool fix = cf.getValueSafe<bool>("setup","fix_mode_amplitude",false);
bool flip = cf.getValueSafe<bool>("setup","flip_mode_amplitude",false);
bool fourier_splicing = cf.getValueSafe<bool>("setup","fourier_splicing",true);
if( fix && levelmin != levelmax ){
LOGWARN("You have chosen mode fixing for a zoom. This is not well tested,\n please proceed at your own risk...");
@ -448,8 +455,6 @@ void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type t
convolution::kernel_creator *the_kernel_creator;
if (kspaceTF)
{
std::cout << " - Using k-space transfer function kernel.\n";
LOGUSER("Using k-space transfer function kernel.");
@ -458,17 +463,6 @@ void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type t
#else
the_kernel_creator = convolution::get_kernel_map()["tf_kernel_k_double"];
#endif
}
else
{
std::cout << " - Using real-space transfer function kernel.\n";
LOGUSER("Using real-space transfer function kernel.");
#ifdef SINGLE_PRECISION
the_kernel_creator = convolution::get_kernel_map()["tf_kernel_real_float"];
#else
the_kernel_creator = convolution::get_kernel_map()["tf_kernel_real_double"];
#endif
}
convolution::kernel *the_tf_kernel = the_kernel_creator->create(cf, ptf, refh, type);
@ -501,21 +495,13 @@ void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type t
refh.size(levelmin + i, 1), refh.size(levelmin + i, 2));
if( refh.get_margin() > 0 ){
fine = new PaddedDensitySubGrid<real_t>(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 = new PaddedDensitySubGrid<real_t>( 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),
refh.get_margin(), refh.get_margin(), refh.get_margin() );
LOGUSER(" margin = %d",refh.get_margin());
}else{
fine = new PaddedDensitySubGrid<real_t>(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 = new PaddedDensitySubGrid<real_t>( 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(" margin = %d",refh.size(levelmin + i, 0)/2);
}
/////////////////////////////////////////////////////////////////////////
@ -526,10 +512,12 @@ void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type t
convolution::perform<real_t>(the_tf_kernel->fetch_kernel(levelmin + i, true),
reinterpret_cast<void *>(fine->get_data_ptr()), shift, fix, flip);
if( fourier_splicing ){
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),
@ -563,6 +551,9 @@ void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type t
std::cout << " - Density calculation took " << tend - tstart << "s." << std::endl;
#endif
if( !fourier_splicing ){
coarsen_density(refh,delta,false);
}
LOGUSER("Finished computing the density field in %fs", tend - tstart);
}

284
src/densities.hh
View file

@ -15,298 +15,20 @@
#include "general.hh"
#include "config_file.hh"
// #include "density_grid.hh"
#include "random.hh"
#include "cosmology.hh"
#include "transfer_function.hh"
#include "general.hh"
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, noise_generator &rand, grid_hierarchy &delta, bool smooth, bool shift);
void GenerateDensityUnigrid(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, noise_generator &rand, grid_hierarchy &delta, bool smooth, bool shift);
void normalize_density(grid_hierarchy &delta);
/*!
* @class DensityGrid
* @brief provides infrastructure for computing the initial density field
*
* This class provides access and data management member functions that
* are used when computing the initial density field by convolution with
* transfer functions.
*/
template <typename real_t>
class DensityGrid
{
public:
size_t nx_; //!< number of grid cells in x-direction
size_t ny_; //!< number of grid cells in y-direction
size_t nz_; //!< number of grid cells in z-direction
size_t nzp_; //!< number of cells in memory (z-dir), used for Nyquist padding
size_t nv_[3];
int ox_; //!< offset of grid in x-direction
int oy_; //!< offset of grid in y-direction
int oz_; //!< offset of grid in z-direction
size_t ov_[3];
//! the actual data container in the form of a 1D array
std::vector<real_t> data_;
//! constructor
/*! constructs an instance given the dimensions of the density field
* @param nx the number of cells in x
* @param ny the number of cells in y
* @param nz the number of cells in z
*/
DensityGrid(unsigned nx, unsigned ny, unsigned nz)
: nx_(nx), ny_(ny), nz_(nz), nzp_(2 * (nz_ / 2 + 1)), ox_(0), oy_(0), oz_(0)
{
data_.assign((size_t)nx_ * (size_t)ny_ * (size_t)nzp_, 0.0);
nv_[0] = nx_;
nv_[1] = ny_;
nv_[2] = nz_;
ov_[0] = ox_;
ov_[1] = oy_;
ov_[2] = oz_;
}
DensityGrid(unsigned nx, unsigned ny, unsigned nz, int ox, int oy, int oz)
: nx_(nx), ny_(ny), nz_(nz), nzp_(2 * (nz_ / 2 + 1)), ox_(ox), oy_(oy), oz_(oz)
{
data_.assign((size_t)nx_ * (size_t)ny_ * (size_t)nzp_, 0.0);
nv_[0] = nx_;
nv_[1] = ny_;
nv_[2] = nz_;
ov_[0] = ox_;
ov_[1] = oy_;
ov_[2] = oz_;
}
//! copy constructor
explicit DensityGrid(const DensityGrid<real_t> &g)
: nx_(g.nx_), ny_(g.ny_), nz_(g.nz_), nzp_(g.nzp_),
ox_(g.ox_), oy_(g.oy_), oz_(g.oz_)
{
data_ = g.data_;
nv_[0] = nx_;
nv_[1] = ny_;
nv_[2] = nz_;
ov_[0] = ox_;
ov_[1] = oy_;
ov_[2] = oz_;
}
//! destructor
~DensityGrid()
{
}
//! clears the density object
/*! sets all dimensions to zero and frees the memory
*/
void clear(void)
{
nx_ = ny_ = nz_ = nzp_ = 0;
ox_ = oy_ = oz_ = 0;
nv_[0] = nv_[1] = nv_[2] = 0;
ov_[0] = ov_[1] = ov_[2] = 0;
data_.clear();
std::vector<real_t>().swap(data_);
}
//! query the 3D array sizes of the density object
/*! returns the size of the 3D density object along a specified dimension
* @param i the dimension for which size is to be returned
* @returns array size along dimension i
*/
size_t size(int i)
{
return nv_[i];
}
int offset(int i)
{
return ov_[i];
}
//! zeroes the density object
/*! sets all values to 0.0
*/
void zero(void)
{
data_.assign(data_.size(), 0.0);
}
//! assigns the contents of another DensityGrid to this
DensityGrid &operator=(const DensityGrid<real_t> &g)
{
nx_ = g.nx_;
ny_ = g.ny_;
nz_ = g.nz_;
nzp_ = g.nzp_;
ox_ = g.ox_;
oy_ = g.oy_;
oz_ = g.oz_;
data_ = g.data_;
return *this;
}
//! 3D index based data access operator
inline real_t &operator()(size_t i, size_t j, size_t k)
{
return data_[((size_t)i * ny_ + (size_t)j) * nzp_ + (size_t)k];
}
//! 3D index based const data access operator
inline const real_t &operator()(size_t i, size_t j, size_t k) const
{
return data_[((size_t)i * ny_ + (size_t)j) * nzp_ + (size_t)k];
}
//! recover the pointer to the 1D data array
inline real_t *get_data_ptr(void)
{
return &data_[0];
}
//! fills the density field with random number values
/*! given a pointer to a random_numbers object, fills the field with random values
* @param prc pointer to a random_numbers object
* @param variance the variance of the random numbers (the values returned by prc are multiplied by this)
* @param i0 x-offset (shift) in cells of the density field with respect to the random number field
* @param j0 y-offset (shift) in cells of the density field with respect to the random number field
* @param k0 z-offset (shift) in cells of the density field with respect to the random number field
* @param setzero boolean, if true, the global mean will be subtracted
*/
void fill_rand(/*const*/ random_numbers<real_t> *prc, real_t variance, int i0, int j0, int k0, bool setzero = false)
{
long double sum = 0.0;
#pragma omp parallel for reduction(+ \
: sum)
for (int i = 0; i < nx_; ++i)
for (int j = 0; j < ny_; ++j)
for (int k = 0; k < nz_; ++k)
{
(*this)(i, j, k) = (*prc)(i0 + i, j0 + j, k0 + k) * variance;
sum += (*this)(i, j, k);
}
sum /= nx_ * ny_ * nz_;
if (setzero)
{
#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)
(*this)(i, j, k) -= sum;
}
}
//! copies the data from another field with 3D index-based access operator
template <class array3>
void copy(array3 &v)
{
#pragma omp parallel for
for (int ix = 0; ix < (int)nx_; ++ix)
for (int iy = 0; iy < (int)ny_; ++iy)
for (int iz = 0; iz < (int)nz_; ++iz)
v(ix, iy, iz) = (*this)(ix, iy, iz);
}
//! adds the data from another field with 3D index-based access operator
template <class array3>
void copy_add(array3 &v)
{
#pragma omp parallel for
for (int ix = 0; ix < (int)nx_; ++ix)
for (int iy = 0; iy < (int)ny_; ++iy)
for (int iz = 0; iz < (int)nz_; ++iz)
v(ix, iy, iz) += (*this)(ix, iy, iz);
}
};
template <typename real_t>
class PaddedDensitySubGrid : public DensityGrid<real_t>
{
public:
using DensityGrid<real_t>::nx_;
using DensityGrid<real_t>::ny_;
using DensityGrid<real_t>::nz_;
using DensityGrid<real_t>::ox_;
using DensityGrid<real_t>::oy_;
using DensityGrid<real_t>::oz_;
using DensityGrid<real_t>::data_;
std::array<size_t,3> pad_;
using DensityGrid<real_t>::fill_rand;
using DensityGrid<real_t>::get_data_ptr;
public:
PaddedDensitySubGrid(int ox, int oy, int oz, unsigned nx, unsigned ny, unsigned nz )
: DensityGrid<real_t>(nx*2, ny*2, nz*2, ox, oy, oz),
pad_{{ nx / 2, ny / 2, nz / 2 }}
{ }
PaddedDensitySubGrid(int ox, int oy, int oz, unsigned nx, unsigned ny, unsigned nz, unsigned padx, unsigned pady, unsigned padz )
: DensityGrid<real_t>(nx + 2 * padx, ny + 2 * pady, nz + 2 * padz, ox, oy, oz),
pad_{{ padx, pady, padz }}
{ }
PaddedDensitySubGrid(const PaddedDensitySubGrid<real_t> &o)
: DensityGrid<real_t>(o)
{ }
size_t margin(int i) const
{
return pad_[i];
}
template <class array3>
void copy_unpad(array3 &v)
{
#pragma omp parallel for
for (size_t ix = pad_[0]; ix < nx_-pad_[0]; ++ix){
const size_t ixu = ix - pad_[0];
for (size_t iy = pad_[1], iyu = 0; iy < ny_-pad_[1]; ++iy, ++iyu)
for (size_t iz = pad_[2], izu = 0; iz < nz_-pad_[2]; ++iz, ++izu)
v(ixu, iyu, izu) = (*this)(ix, iy, iz);
}
}
template <class array3>
void copy_add_unpad(array3 &v)
{
#pragma omp parallel for
for (size_t ix = pad_[0]; ix < nx_-pad_[0]; ++ix){
const size_t ixu = ix - pad_[0];
for (size_t iy = pad_[1], iyu = 0; iy < ny_-pad_[1]; ++iy, ++iyu)
for (size_t iz = pad_[2], izu = 0; iz < nz_-pad_[2]; ++iz, ++izu)
v(ixu, iyu, izu) += (*this)(ix, iy, iz);
}
}
template <class array3>
void copy_subtract_unpad(array3 &v)
{
#pragma omp parallel for
for (size_t ix = pad_[0]; ix < nx_-pad_[0]; ++ix){
const size_t ixu = ix - pad_[0];
for (size_t iy = pad_[1], iyu = 0; iy < ny_-pad_[1]; ++iy, ++iyu)
for (size_t iz = pad_[2], izu = 0; iz < nz_-pad_[2]; ++iz, ++izu)
v(ixu, iyu, izu) -= (*this)(ix, iy, iz);
}
}
};
void coarsen_density(const refinement_hierarchy &rh, GridHierarchy<real_t> &u, bool kspace);
#endif

297
src/density_grid.hh Normal file
View file

@ -0,0 +1,297 @@
#pragma once
#include <vector>
#include <array>
/*!
* @class DensityGrid
* @brief provides infrastructure for computing the initial density field
*
* This class provides access and data management member functions that
* are used when computing the initial density field by convolution with
* transfer functions.
*/
template <typename real_t>
class DensityGrid
{
public:
size_t nx_; //!< number of grid cells in x-direction
size_t ny_; //!< number of grid cells in y-direction
size_t nz_; //!< number of grid cells in z-direction
size_t nzp_; //!< number of cells in memory (z-dir), used for Nyquist padding
std::array<size_t,3> nv_;
int ox_; //!< offset of grid in x-direction
int oy_; //!< offset of grid in y-direction
int oz_; //!< offset of grid in z-direction
std::array<int,3> ov_;
//! the actual data container in the form of a 1D array
std::vector<real_t> data_;
//! constructor
/*! constructs an instance given the dimensions of the density field
* @param nx the number of cells in x
* @param ny the number of cells in y
* @param nz the number of cells in z
*/
DensityGrid(unsigned nx, unsigned ny, unsigned nz)
: nx_(nx), ny_(ny), nz_(nz), nzp_(2 * (nz_ / 2 + 1)), ox_(0), oy_(0), oz_(0)
{
data_.assign((size_t)nx_ * (size_t)ny_ * (size_t)nzp_, 0.0);
nv_ = {nx_,ny_,nz_};
ov_ = {ox_,oy_,oz_};
// nv_[0] = nx_;
// nv_[1] = ny_;
// nv_[2] = nz_;
// ov_[0] = ox_;
// ov_[1] = oy_;
// ov_[2] = oz_;
}
DensityGrid(unsigned nx, unsigned ny, unsigned nz, int ox, int oy, int oz)
: nx_(nx), ny_(ny), nz_(nz), nzp_(2 * (nz_ / 2 + 1)), ox_(ox), oy_(oy), oz_(oz)
{
data_.assign((size_t)nx_ * (size_t)ny_ * (size_t)nzp_, 0.0);
nv_ = {nx_,ny_,nz_};
ov_ = {ox_,oy_,oz_};
// nv_[0] = nx_;
// nv_[1] = ny_;
// nv_[2] = nz_;
// ov_[0] = ox_;
// ov_[1] = oy_;
// ov_[2] = oz_;
}
//! copy constructor
explicit DensityGrid(const DensityGrid<real_t> &g)
: nx_(g.nx_), ny_(g.ny_), nz_(g.nz_), nzp_(g.nzp_),
ox_(g.ox_), oy_(g.oy_), oz_(g.oz_)
{
data_ = g.data_;
nv_ = {nx_,ny_,nz_};
ov_ = {ox_,oy_,oz_};
// nv_[0] = nx_;
// nv_[1] = ny_;
// nv_[2] = nz_;
// ov_[0] = ox_;
// ov_[1] = oy_;
// ov_[2] = oz_;
}
//! destructor
~DensityGrid()
{
}
//! clears the density object
/*! sets all dimensions to zero and frees the memory
*/
void clear(void)
{
nx_ = ny_ = nz_ = nzp_ = 0;
ox_ = oy_ = oz_ = 0;
nv_[0] = nv_[1] = nv_[2] = 0;
ov_[0] = ov_[1] = ov_[2] = 0;
data_.clear();
std::vector<real_t>().swap(data_);
}
//! query the 3D array sizes of the density object
/*! returns the size of the 3D density object along a specified dimension
* @param i the dimension for which size is to be returned
* @returns array size along dimension i
*/
size_t size(int i)
{
return nv_[i];
}
int offset(int i)
{
return ov_[i];
}
//! zeroes the density object
/*! sets all values to 0.0
*/
void zero(void)
{
data_.assign(data_.size(), 0.0);
}
//! assigns the contents of another DensityGrid to this
DensityGrid &operator=(const DensityGrid<real_t> &g)
{
nx_ = g.nx_;
ny_ = g.ny_;
nz_ = g.nz_;
nzp_ = g.nzp_;
ox_ = g.ox_;
oy_ = g.oy_;
oz_ = g.oz_;
data_ = g.data_;
return *this;
}
//! 3D index based data access operator
inline real_t &operator()(size_t i, size_t j, size_t k)
{
return data_[((size_t)i * ny_ + (size_t)j) * nzp_ + (size_t)k];
}
//! 3D index based const data access operator
inline const real_t &operator()(size_t i, size_t j, size_t k) const
{
return data_[((size_t)i * ny_ + (size_t)j) * nzp_ + (size_t)k];
}
//! recover the pointer to the 1D data array
inline real_t *get_data_ptr(void)
{
return &data_[0];
}
#if 0
//! fills the density field with random number values
/*! given a pointer to a random_numbers object, fills the field with random values
* @param prc pointer to a random_numbers object
* @param variance the variance of the random numbers (the values returned by prc are multiplied by this)
* @param i0 x-offset (shift) in cells of the density field with respect to the random number field
* @param j0 y-offset (shift) in cells of the density field with respect to the random number field
* @param k0 z-offset (shift) in cells of the density field with respect to the random number field
* @param setzero boolean, if true, the global mean will be subtracted
*/
void fill_rand(/*const*/ random_numbers<real_t> *prc, real_t variance, int i0, int j0, int k0, bool setzero = false)
{
long double sum = 0.0;
#pragma omp parallel for reduction(+ \
: sum)
for (int i = 0; i < nx_; ++i)
for (int j = 0; j < ny_; ++j)
for (int k = 0; k < nz_; ++k)
{
(*this)(i, j, k) = (*prc)(i0 + i, j0 + j, k0 + k) * variance;
sum += (*this)(i, j, k);
}
sum /= nx_ * ny_ * nz_;
if (setzero)
{
#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)
(*this)(i, j, k) -= sum;
}
}
#endif
//! copies the data from another field with 3D index-based access operator
template <class array3>
void copy(array3 &v)
{
#pragma omp parallel for
for (int ix = 0; ix < (int)nx_; ++ix)
for (int iy = 0; iy < (int)ny_; ++iy)
for (int iz = 0; iz < (int)nz_; ++iz)
v(ix, iy, iz) = (*this)(ix, iy, iz);
}
//! adds the data from another field with 3D index-based access operator
template <class array3>
void copy_add(array3 &v)
{
#pragma omp parallel for
for (int ix = 0; ix < (int)nx_; ++ix)
for (int iy = 0; iy < (int)ny_; ++iy)
for (int iz = 0; iz < (int)nz_; ++iz)
v(ix, iy, iz) += (*this)(ix, iy, iz);
}
};
template <typename real_t>
class PaddedDensitySubGrid : public DensityGrid<real_t>
{
public:
using DensityGrid<real_t>::nx_;
using DensityGrid<real_t>::ny_;
using DensityGrid<real_t>::nz_;
using DensityGrid<real_t>::ox_;
using DensityGrid<real_t>::oy_;
using DensityGrid<real_t>::oz_;
using DensityGrid<real_t>::data_;
std::array<size_t, 3> pad_;
// using DensityGrid<real_t>::fill_rand;
using DensityGrid<real_t>::get_data_ptr;
public:
PaddedDensitySubGrid(int ox, int oy, int oz, unsigned nx, unsigned ny, unsigned nz)
: DensityGrid<real_t>(nx * 2, ny * 2, nz * 2, ox, oy, oz),
pad_{{nx / 2, ny / 2, nz / 2}}
{
}
PaddedDensitySubGrid(int ox, int oy, int oz, unsigned nx, unsigned ny, unsigned nz, unsigned padx, unsigned pady, unsigned padz)
: DensityGrid<real_t>(nx + 2 * padx, ny + 2 * pady, nz + 2 * padz, ox, oy, oz),
pad_{{padx, pady, padz}}
{
}
PaddedDensitySubGrid(const PaddedDensitySubGrid<real_t> &o)
: DensityGrid<real_t>(o)
{
}
size_t margin(int i) const
{
return pad_[i];
}
template <class array3>
void copy_unpad(array3 &v)
{
#pragma omp parallel for
for (size_t ix = pad_[0]; ix < nx_ - pad_[0]; ++ix)
{
const size_t ixu = ix - pad_[0];
for (size_t iy = pad_[1], iyu = 0; iy < ny_ - pad_[1]; ++iy, ++iyu)
for (size_t iz = pad_[2], izu = 0; iz < nz_ - pad_[2]; ++iz, ++izu)
v(ixu, iyu, izu) = (*this)(ix, iy, iz);
}
}
template <class array3>
void copy_add_unpad(array3 &v)
{
#pragma omp parallel for
for (size_t ix = pad_[0]; ix < nx_ - pad_[0]; ++ix)
{
const size_t ixu = ix - pad_[0];
for (size_t iy = pad_[1], iyu = 0; iy < ny_ - pad_[1]; ++iy, ++iyu)
for (size_t iz = pad_[2], izu = 0; iz < nz_ - pad_[2]; ++iz, ++izu)
v(ixu, iyu, izu) += (*this)(ix, iy, iz);
}
}
template <class array3>
void copy_subtract_unpad(array3 &v)
{
#pragma omp parallel for
for (size_t ix = pad_[0]; ix < nx_ - pad_[0]; ++ix)
{
const size_t ixu = ix - pad_[0];
for (size_t iy = pad_[1], iyu = 0; iy < ny_ - pad_[1]; ++iy, ++iyu)
for (size_t iz = pad_[2], izu = 0; iz < nz_ - pad_[2]; ++iz, ++izu)
v(ixu, iyu, izu) -= (*this)(ix, iy, iz);
}
}
};

204
src/fft_operators.hh
View file

@ -1,204 +0,0 @@
#ifndef __FFT_OPERATORS_HH
#define __FFT_OPERATORS_HH
struct fft_interp{
template< typename m1, typename m2 >
void interpolate( m1& V, m2& v, bool fourier_splice = false ) const
{
int oxc = V.offset(0), oyc = V.offset(1), ozc = V.offset(2);
int oxf = v.offset(0), oyf = v.offset(1), ozf = v.offset(2);
size_t nxf = v.size(0), nyf = v.size(1), nzf = v.size(2), nzfp = nzf+2;
// cut out piece of coarse grid that overlaps the fine:
assert( nxf%2==0 && nyf%2==0 && nzf%2==0 );
size_t nxc = nxf/2, nyc = nyf/2, nzc = nzf/2, nzcp = nzf/2+2;
fftw_real *rcoarse = new fftw_real[ nxc * nyc * nzcp ];
fftw_complex *ccoarse = reinterpret_cast<fftw_complex*> (rcoarse);
fftw_real *rfine = new fftw_real[ nxf * nyf * nzfp];
fftw_complex *cfine = reinterpret_cast<fftw_complex*> (rfine);
#pragma omp parallel for
for( int i=0; i<(int)nxc; ++i )
for( int j=0; j<(int)nyc; ++j )
for( int k=0; k<(int)nzc; ++k )
{
size_t q = ((size_t)i*nyc+(size_t)j)*nzcp+(size_t)k;
rcoarse[q] = V( oxf+i, oyf+j, ozf+k );
}
if( fourier_splice )
{
#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);
}
}
else
{
#pragma omp parallel for
for( size_t i=0; i<nxf*nyf*nzfp; ++i )
rfine[i] = 0.0;
}
#ifdef FFTW3
#ifdef SINGLE_PRECISION
fftwf_plan
pc = fftwf_plan_dft_r2c_3d( nxc, nyc, nzc, rcoarse, ccoarse, FFTW_ESTIMATE),
pf = fftwf_plan_dft_r2c_3d( nxf, nyf, nzf, rfine, cfine, FFTW_ESTIMATE),
ipf = fftwf_plan_dft_c2r_3d( nxf, nyf, nzf, cfine, rfine, FFTW_ESTIMATE);
fftwf_execute( pc );
if( fourier_splice )
fftwf_execute( pf );
#else
fftw_plan
pc = fftw_plan_dft_r2c_3d( nxc, nyc, nzc, rcoarse, ccoarse, FFTW_ESTIMATE),
pf = fftw_plan_dft_r2c_3d( nxf, nyf, nzf, rfine, cfine, FFTW_ESTIMATE),
ipf = fftw_plan_dft_c2r_3d( nxf, nyf, nzf, cfine, rfine, FFTW_ESTIMATE);
fftw_execute( pc );
if( fourier_splice )
fftwf_execute( pf );
#endif
#else
rfftwnd_plan
pc = rfftw3d_create_plan( nxc, nyc, nzc, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE|FFTW_IN_PLACE),
pf = rfftw3d_create_plan( nxf, nyf, nzf, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE|FFTW_IN_PLACE),
ipf = rfftw3d_create_plan( nxf, nyf, nzf, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE|FFTW_IN_PLACE);
#ifndef SINGLETHREAD_FFTW
rfftwnd_threads_one_real_to_complex( omp_get_max_threads(), pc, rcoarse, NULL );
if( fourier_splice )
rfftwnd_threads_one_real_to_complex( omp_get_max_threads(), pf, rfine, NULL );
#else
rfftwnd_one_real_to_complex( pc, rcoarse, NULL );
if( fourier_splice )
rfftwnd_one_real_to_complex( pf, rfine, NULL );
#endif
#endif
/*************************************************/
//.. perform actual interpolation
double fftnorm = 1.0/((double)nxf*(double)nyf*(double)nzf);
double sqrt8 = sqrt(8.0);
// 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);
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)nyf+(size_t)jj)*(nzf/2+1)+(size_t)kk;
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( k==0 & (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( k==0 && (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)nyf+(size_t)jj)*(nzf/2+1)+(size_t)kk;
RE(cfine[qf]) = sqrt8*RE(ccoarse[qc]);
IM(cfine[qf]) = sqrt8*IM(ccoarse[qc]);
}
delete[] rcoarse;
/*************************************************/
#ifdef FFTW3
#ifdef SINGLE_PRECISION
fftwf_execute( ipf );
fftwf_destroy_plan(pf);
fftwf_destroy_plan(pc);
fftwf_destroy_plan(ipf);
#else
fftw_execute( ipf );
fftw_destroy_plan(pf);
fftw_destroy_plan(pc);
fftw_destroy_plan(ipf);
#endif
#else
#ifndef SINGLETHREAD_FFTW
rfftwnd_threads_one_complex_to_real( omp_get_max_threads(), ipf, cfine, NULL );
#else
rfftwnd_one_complex_to_real( ipf, cfine, NULL );
#endif
fftwnd_destroy_plan(pf);
fftwnd_destroy_plan(pc);
fftwnd_destroy_plan(ipf);
#endif
// copy back and normalize
#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;
v(i,j,k) = rfine[q] * fftnorm;
}
delete[] rcoarse;
delete[] rfine;
}
};
#endif //__FFT_OPERATORS_HH

52
src/main.cc
View file

@ -43,7 +43,7 @@ extern "C"
#include "transfer_function.hh"
#define THE_CODE_NAME "music!"
#define THE_CODE_VERSION "1.53"
#define THE_CODE_VERSION "2.0a"
namespace music
{
@ -249,7 +249,7 @@ double compute_finest_sigma(grid_hierarchy &u)
return sqrt(sum2 - sum * sum);
}
double compute_finest_max(grid_hierarchy &u)
double compute_finest_absmax(grid_hierarchy &u)
{
double valmax = 0.0;
#pragma omp parallel for reduction(max:valmax)
@ -257,8 +257,8 @@ double compute_finest_max(grid_hierarchy &u)
for (int iy = 0; iy < (int)(*u.get_grid(u.levelmax())).size(1); ++iy)
for (int iz = 0; iz < (int)(*u.get_grid(u.levelmax())).size(2); ++iz)
{
if (fabs((*u.get_grid(u.levelmax()))(ix, iy, iz)) > fabs(valmax))
valmax = (*u.get_grid(u.levelmax()))(ix, iy, iz);
if (std::fabs((*u.get_grid(u.levelmax()))(ix, iy, iz)) > valmax)
valmax = std::fabs((*u.get_grid(u.levelmax()))(ix, iy, iz));
}
return valmax;
@ -299,7 +299,6 @@ void add_constant_value( grid_hierarchy &u, const double val )
/*****************************************************************************************************/
region_generator_plugin *the_region_generator;
RNG_plugin *the_random_number_generator;
int main(int argc, const char *argv[])
{
@ -422,7 +421,8 @@ int main(int argc, const char *argv[])
bool
do_baryons = cf.getValue<bool>("setup", "baryons"),
do_2LPT = cf.getValueSafe<bool>("setup", "use_2LPT", false),
do_LLA = cf.getValueSafe<bool>("setup", "use_LLA", false);
do_LLA = cf.getValueSafe<bool>("setup", "use_LLA", false),
do_counter_mode = cf.getValueSafe<bool>("setup", "zero_zoom_velocity", false);
transfer_function_plugin *the_transfer_function_plugin = select_transfer_function_plugin(cf);
@ -462,7 +462,6 @@ int main(int argc, const char *argv[])
the_region_generator = select_region_generator_plugin(cf);
the_random_number_generator = select_RNG_plugin(cf);
//------------------------------------------------------------------------------
//... determine run parameters
//------------------------------------------------------------------------------
@ -472,6 +471,12 @@ int main(int argc, const char *argv[])
<< " distinct amplitudes for baryon and DM fields!\n"
<< " Perturbation amplitudes will be identical!" << std::endl;
//------------------------------------------------------------------------------
//... start up the random number generator plugin
//... see if we need to set some grid building constraints
noise_generator rand( cf, the_transfer_function_plugin );
//------------------------------------------------------------------------------
//... determine the refinement hierarchy
//------------------------------------------------------------------------------
@ -505,7 +510,8 @@ int main(int argc, const char *argv[])
std::cout << " GENERATING WHITE NOISE\n";
std::cout << "-------------------------------------------------------------\n";
LOGUSER("Computing white noise...");
rand_gen rand(cf, rh_TF, the_transfer_function_plugin);
// rand_gen rand(cf, rh_TF, the_transfer_function_plugin);
rand.initialize_for_grid_structure( rh_TF );
//------------------------------------------------------------------------------
//... initialize the Poisson solver
@ -608,13 +614,13 @@ int main(int argc, const char *argv[])
else
//... displacement
the_poisson_solver->gradient(icoord, u, data_forIO);
double dispmax = compute_finest_max(data_forIO);
double dispmax = compute_finest_absmax(data_forIO);
LOGINFO("max. %c-displacement of HR particles is %f [mean dx]", 'x' + icoord, dispmax * (double)(1ll << data_forIO.levelmax()));
coarsen_density(rh_Poisson, data_forIO, false);
//... compute counter-mode to minimize advection errors
counter_mode_amp[icoord] = compute_finest_mean(data_forIO);
add_constant_value( data_forIO, -counter_mode_amp[icoord] );
if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord] );
LOGUSER("Writing CDM displacements");
the_output_plugin->write_dm_position(icoord, data_forIO);
@ -745,13 +751,13 @@ int main(int argc, const char *argv[])
LOGINFO("mean of %c-velocity of high-res particles is %f", 'x' + icoord, meanv);
LOGUSER("mean of %c-velocity of high-res particles is %f", 'x' + icoord, meanv);
double maxv = compute_finest_max(data_forIO);
double maxv = compute_finest_absmax(data_forIO);
LOGINFO("max of abs of %c-velocity of high-res particles is %f", 'x' + icoord, maxv);
coarsen_density(rh_Poisson, data_forIO, false);
// add counter velocity-mode
add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact );
if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact );
LOGUSER("Writing CDM velocities");
the_output_plugin->write_dm_velocity(icoord, data_forIO);
@ -815,13 +821,13 @@ int main(int argc, const char *argv[])
LOGINFO("mean of %c-velocity of high-res particles is %f", 'x' + icoord, meanv);
LOGUSER("mean of %c-velocity of high-res particles is %f", 'x' + icoord, meanv);
double maxv = compute_finest_max(data_forIO);
double maxv = compute_finest_absmax(data_forIO);
LOGINFO("max of abs of %c-velocity of high-res particles is %f", 'x' + icoord, maxv);
coarsen_density(rh_Poisson, data_forIO, false);
// add counter velocity mode
add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact );
if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact );
LOGUSER("Writing CDM velocities");
the_output_plugin->write_dm_velocity(icoord, data_forIO);
@ -874,13 +880,13 @@ int main(int argc, const char *argv[])
LOGINFO("mean of %c-velocity of high-res baryons is %f", 'x' + icoord, meanv);
LOGUSER("mean of %c-velocity of high-res baryons is %f", 'x' + icoord, meanv);
double maxv = compute_finest_max(data_forIO);
double maxv = compute_finest_absmax(data_forIO);
LOGINFO("max of abs of %c-velocity of high-res baryons is %f", 'x' + icoord, maxv);
coarsen_density(rh_Poisson, data_forIO, false);
// add counter velocity mode
add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact );
if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact );
LOGUSER("Writing baryon velocities");
the_output_plugin->write_gas_velocity(icoord, data_forIO);
@ -993,7 +999,7 @@ int main(int argc, const char *argv[])
LOGINFO("mean of %c-velocity of high-res particles is %f", 'x' + icoord, meanv);
LOGUSER("mean of %c-velocity of high-res particles is %f", 'x' + icoord, meanv);
double maxv = compute_finest_max(data_forIO);
double maxv = compute_finest_absmax(data_forIO);
LOGINFO("max of abs of %c-velocity of high-res particles is %f", 'x' + icoord, maxv);
std::cerr << " - velocity component " << icoord << " : sigma = " << sigv << std::endl;
@ -1003,7 +1009,7 @@ int main(int argc, const char *argv[])
//... compute counter-mode to minimize advection errors
counter_mode_amp[icoord] = compute_finest_mean(data_forIO);
add_constant_value( data_forIO, -counter_mode_amp[icoord] );
if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord] );
LOGUSER("Writing CDM velocities");
the_output_plugin->write_dm_velocity(icoord, data_forIO);
@ -1091,7 +1097,7 @@ int main(int argc, const char *argv[])
LOGINFO("mean of %c-velocity of high-res baryons is %f", 'x' + icoord, meanv);
LOGUSER("mean of %c-velocity of high-res baryons is %f", 'x' + icoord, meanv);
double maxv = compute_finest_max(data_forIO);
double maxv = compute_finest_absmax(data_forIO);
LOGINFO("max of abs of %c-velocity of high-res baryons is %f", 'x' + icoord, maxv);
std::cerr << " - velocity component " << icoord << " : sigma = " << sigv << std::endl;
@ -1100,7 +1106,7 @@ int main(int argc, const char *argv[])
coarsen_density(rh_Poisson, data_forIO, false);
// add counter velocity mode
add_constant_value( data_forIO, -counter_mode_amp[icoord] );
if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord] );
LOGUSER("Writing baryon velocities");
the_output_plugin->write_gas_velocity(icoord, data_forIO);
@ -1199,13 +1205,13 @@ int main(int argc, const char *argv[])
else
the_poisson_solver->gradient(icoord, u1, data_forIO);
double dispmax = compute_finest_max(data_forIO);
double dispmax = compute_finest_absmax(data_forIO);
LOGINFO("max. %c-displacement of HR particles is %f [mean dx]", 'x' + icoord, dispmax * (double)(1ll << data_forIO.levelmax()));
coarsen_density(rh_Poisson, data_forIO, false);
// add counter mode
add_constant_value( data_forIO, -counter_mode_amp[icoord]/cosmo.vfact );
if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]/cosmo.vfact );
LOGUSER("Writing CDM displacements");
the_output_plugin->write_dm_position(icoord, data_forIO);
@ -1322,7 +1328,7 @@ int main(int argc, const char *argv[])
coarsen_density(rh_Poisson, data_forIO, false);
// add counter mode
add_constant_value( data_forIO, -counter_mode_amp[icoord]/cosmo.vfact );
if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]/cosmo.vfact );
LOGUSER("Writing baryon displacements");

81
src/mesh.hh
View file

@ -1212,7 +1212,8 @@ class refinement_hierarchy
levelmax_, //!< maximum grid level for all operations
levelmin_tf_, //!< minimum grid level for density calculation
padding_, //!< padding in number of coarse cells between refinement levels
blocking_factor_;
blocking_factor_, //!< blocking factor of grids, necessary fo BoxLib codes such as NyX
gridding_unit_; //!< internal blocking factor of grids, necessary for Panphasia
int margin_; //!< number of cells used for additional padding for convolutions with isolated boundaries (-1 = double padding)
@ -1231,6 +1232,30 @@ class refinement_hierarchy
index3_t xshift_; //!< shift of refinement region in coarse cells (in order to center it in the domain)
double rshift_[3];
//! calculates the greatest common divisor
int gcd(int a, int b) const {
return b == 0 ? a : gcd(b, a % b);
}
// calculates the cell shift in units of levelmin grid cells if there is an additional constraint to be
// congruent with another grid that partitions the same space in multiples of "base_unit"
int get_shift_unit( int base_unit, int levelmin ) const {
/*int Lp = 0;
while( base_unit * (1<<Lp) < 1<<(levelmin+1) ){
++Lp;
}
int U = base_unit * (1<<Lp);
return std::max<int>( 1, (1<<(levelmin+1)) / (2*gcd(U,1<<(levelmin+1) )) );*/
int level_m = 0;
while( base_unit * (1<<level_m) < (1<<levelmin) )
++level_m;
return std::max<int>( 1, (1<<levelmin)/gcd(base_unit * (1<<level_m),(1<<levelmin)) );
}
public:
//! copy constructor
refinement_hierarchy(const refinement_hierarchy &rh)
@ -1256,6 +1281,16 @@ public:
bool bnoshift = cf_.getValueSafe<bool>("setup", "no_shift", false);
bool force_shift = cf_.getValueSafe<bool>("setup", "force_shift", false);
gridding_unit_ = cf.getValueSafe<unsigned>("setup", "gridding_unit", 2);
if (gridding_unit_ != 2 && blocking_factor_==0) {
blocking_factor_ = gridding_unit_; // THIS WILL LIKELY CAUSE PROBLEMS WITH NYX
}else if (gridding_unit_ != 2 && blocking_factor_!=0 && gridding_unit_!=blocking_factor_ ) {
LOGERR("incompatible gridding unit %d and blocking factor specified", gridding_unit_, blocking_factor_ );
throw std::runtime_error("Incompatible gridding unit and blocking factor!");
}
//... call the region generator
if (levelmin_ != levelmax_)
{
@ -1292,9 +1327,18 @@ public:
if ((levelmin_ != levelmax_) && (!bnoshift || force_shift))
{
xshift_[0] = (int)((0.5 - xc[0]) * ncoarse);
xshift_[1] = (int)((0.5 - xc[1]) * ncoarse);
xshift_[2] = (int)((0.5 - xc[2]) * ncoarse);
int random_base_grid_unit = cf.getValueSafe<int>("random","base_unit",1);
int shift_unit = get_shift_unit( random_base_grid_unit, levelmin_ );
if( shift_unit != 1 ){
LOGINFO("volume can only be shifted by multiples of %d coarse cells.",shift_unit);
}
xshift_[0] = (int)((0.5-xc[0]) * (double)ncoarse / shift_unit + 0.5) * shift_unit;//ARJ(int)((0.5 - xc[0]) * ncoarse);
xshift_[1] = (int)((0.5-xc[1]) * (double)ncoarse / shift_unit + 0.5) * shift_unit;//ARJ(int)((0.5 - xc[1]) * ncoarse);
xshift_[2] = (int)((0.5-xc[2]) * (double)ncoarse / shift_unit + 0.5) * shift_unit;//ARJ(int)((0.5 - xc[2]) * ncoarse);
// xshift_[0] = (int)((0.5 - xc[0]) * ncoarse);
// xshift_[1] = (int)((0.5 - xc[1]) * ncoarse);
// xshift_[2] = (int)((0.5 - xc[2]) * ncoarse);
}
else
{
@ -1414,12 +1458,15 @@ public:
else
{
//... require alignment with coarser grid
il -= il % 2;
jl -= jl % 2;
kl -= kl % 2;
ir += ir % 2;
jr += jr % 2;
kr += kr % 2;
LOGINFO("Internal refinement bounding box error: [%d,%d]x[%d,%d]x[%d,%d]", il, ir, jl, jr, kl, kr);
il -= il % gridding_unit_;
jl -= jl % gridding_unit_;
kl -= kl % gridding_unit_;
ir = ((ir%gridding_unit_)!=0)? (ir/gridding_unit_ + 1)*gridding_unit_ : ir;
jr = ((jr%gridding_unit_)!=0)? (jr/gridding_unit_ + 1)*gridding_unit_ : jr;
kr = ((kr%gridding_unit_)!=0)? (kr/gridding_unit_ + 1)*gridding_unit_ : kr;
}
// if doing unigrid, set region to whole box
@ -1511,7 +1558,6 @@ public:
jr = (int)((double)jr / nref + 1.0) * nref;
kr = (int)((double)kr / nref + 1.0) * nref;
}
else if (preserve_dims_)
{
//... require alignment with coarser grid
@ -1528,12 +1574,13 @@ public:
else
{
//... require alignment with coarser grid
il -= il % 2;
jl -= jl % 2;
kl -= kl % 2;
ir += ir % 2;
jr += jr % 2;
kr += kr % 2;
il -= il % gridding_unit_;
jl -= jl % gridding_unit_;
kl -= kl % gridding_unit_;
ir = ((ir%gridding_unit_)!=0)? (ir/gridding_unit_ + 1)*gridding_unit_ : ir;
jr = ((jr%gridding_unit_)!=0)? (jr/gridding_unit_ + 1)*gridding_unit_ : jr;
kr = ((kr%gridding_unit_)!=0)? (kr/gridding_unit_ + 1)*gridding_unit_ : kr;
}
if (il >= ir || jl >= jr || kl >= kr || il < 0 || jl < 0 || kl < 0)

611
src/plugins/random_music.cc
View file

@ -1,20 +1,609 @@
#include "random.hh"
#include "random_music_wnoise_generator.hh"
typedef music_wnoise_generator<real_t> rng;
class RNG_music : public RNG_plugin
{
protected:
std::vector<long> rngseeds_;
std::vector<std::string> rngfnames_;
unsigned ran_cube_size_;
int levelmin_, levelmax_, levelmin_seed_;
bool disk_cached_;
bool restart_;
bool initialized_;
std::vector<std::vector<real_t> *> mem_cache_;
//! checks if the specified string is numeric
bool is_number(const std::string &s);
//! parses the random number parameters in the conf file
void parse_random_parameters(void);
//! computes the white noise fields and keeps them either in memory or on disk
void compute_random_numbers(void);
//! adjusts averages
void correct_avg(int icoarse, int ifine);
//! store the white noise fields in memory or on disk
void store_rnd(int ilevel, rng *prng);
class RNG_music : public RNG_plugin{
public:
explicit RNG_music( config_file& cf )
: RNG_plugin( cf )
{ }
explicit RNG_music(config_file &cf) : RNG_plugin(cf), initialized_(false) {}
~RNG_music() { }
~RNG_music() {}
bool is_multiscale() const
{ return true; }
bool is_multiscale() const { return true; }
void initialize_for_grid_structure(const refinement_hierarchy &refh)
{
prefh_ = &refh;
levelmin_ = prefh_->levelmin();
levelmax_ = prefh_->levelmax();
ran_cube_size_ = pcf_->getValueSafe<unsigned>("random", "cubesize", DEF_RAN_CUBE_SIZE);
disk_cached_ = pcf_->getValueSafe<bool>("random", "disk_cached", true);
restart_ = pcf_->getValueSafe<bool>("random", "restart", false);
pcf_->insertValue("setup","fourier_splicing","true");
mem_cache_.assign(levelmax_ - levelmin_ + 1, (std::vector<real_t> *)NULL);
if (restart_ && !disk_cached_)
{
LOGERR("Cannot restart from mem cached random numbers.");
throw std::runtime_error("Cannot restart from mem cached random numbers.");
}
//... determine seed/white noise file data to be applied
parse_random_parameters();
if (!restart_)
{
//... compute the actual random numbers
compute_random_numbers();
}
initialized_ = true;
}
void fill_grid(int level, DensityGrid<real_t> &R);
};
bool RNG_music::is_number(const std::string &s)
{
for (size_t i = 0; i < s.length(); i++)
if (!std::isdigit(s[i]) && s[i] != '-')
return false;
namespace{
RNG_plugin_creator_concrete< RNG_music > creator("MUSIC");
return true;
}
void RNG_music::parse_random_parameters(void)
{
//... parse random number options
for (int i = 0; i <= 100; ++i)
{
char seedstr[128];
std::string tempstr;
bool noseed = false;
sprintf(seedstr, "seed[%d]", i);
if (pcf_->containsKey("random", seedstr))
tempstr = pcf_->getValue<std::string>("random", seedstr);
else
{
// "-2" means that no seed entry was found for that level
tempstr = std::string("-2");
noseed = true;
}
if (is_number(tempstr))
{
long ltemp;
pcf_->convert(tempstr, ltemp);
rngfnames_.push_back("");
if (noseed) // ltemp < 0 )
//... generate some dummy seed which only depends on the level, negative so we know it's not
//... an actual seed and thus should not be used as a constraint for coarse levels
// rngseeds_.push_back( -abs((unsigned)(ltemp-i)%123+(unsigned)(ltemp+827342523521*i)%123456789) );
rngseeds_.push_back(-abs((long)(ltemp - i) % 123 + (long)(ltemp + 7342523521 * i) % 123456789));
else
{
if (ltemp <= 0)
{
LOGERR("Specified seed [random]/%s needs to be a number >0!", seedstr);
throw std::runtime_error("Seed values need to be >0");
}
rngseeds_.push_back(ltemp);
}
}
else
{
rngfnames_.push_back(tempstr);
rngseeds_.push_back(-1);
LOGINFO("Random numbers for level %3d will be read from file.", i);
}
}
//.. determine for which levels random seeds/random number files are given
levelmin_seed_ = -1;
for (unsigned ilevel = 0; ilevel < rngseeds_.size(); ++ilevel)
{
if (levelmin_seed_ < 0 && (rngfnames_[ilevel].size() > 0 || rngseeds_[ilevel] >= 0))
levelmin_seed_ = ilevel;
}
}
void RNG_music::compute_random_numbers(void)
{
bool rndsign = pcf_->getValueSafe<bool>("random", "grafic_sign", false);
std::vector<rng *> randc(std::max(levelmax_, levelmin_seed_) + 1, (rng *)NULL);
//--- FILL ALL WHITE NOISE ARRAYS FOR WHICH WE NEED THE FULL FIELD ---//
//... seeds are given for a level coarser than levelmin
if (levelmin_seed_ < levelmin_)
{
if (rngfnames_[levelmin_seed_].size() > 0)
randc[levelmin_seed_] = new rng(1 << levelmin_seed_, rngfnames_[levelmin_seed_], rndsign);
else
randc[levelmin_seed_] = new rng(1 << levelmin_seed_, ran_cube_size_, rngseeds_[levelmin_seed_]);
for (int i = levelmin_seed_ + 1; i <= levelmin_; ++i)
{
// #warning add possibility to read noise from file also here!
if (rngfnames_[i].size() > 0)
LOGINFO("Warning: Cannot use filenames for higher levels currently! Ignoring!");
randc[i] = new rng(*randc[i - 1], ran_cube_size_, rngseeds_[i]);
delete randc[i - 1];
randc[i - 1] = NULL;
}
}
//... seeds are given for a level finer than levelmin, obtain by averaging
if (levelmin_seed_ > levelmin_)
{
if (rngfnames_[levelmin_seed_].size() > 0)
randc[levelmin_seed_] = new rng(1 << levelmin_seed_, rngfnames_[levelmin_seed_], rndsign);
else
randc[levelmin_seed_] =
new rng(1 << levelmin_seed_, ran_cube_size_, rngseeds_[levelmin_seed_]); //, x0, lx );
for (int ilevel = levelmin_seed_ - 1; ilevel >= (int)levelmin_; --ilevel)
{
if (rngseeds_[ilevel - levelmin_] > 0)
LOGINFO("Warning: random seed for level %d will be ignored.\n"
" consistency requires that it is obtained by restriction from level %d",
ilevel, levelmin_seed_);
randc[ilevel] = new rng(*randc[ilevel + 1]);
if (ilevel + 1 > levelmax_)
{
delete randc[ilevel + 1];
randc[ilevel + 1] = NULL;
}
}
}
//--- GENERATE AND STORE ALL LEVELS, INCLUDING REFINEMENTS ---//
//... levelmin
if (randc[levelmin_] == NULL)
{
if (rngfnames_[levelmin_].size() > 0)
randc[levelmin_] = new rng(1 << levelmin_, rngfnames_[levelmin_], rndsign);
else
randc[levelmin_] = new rng(1 << levelmin_, ran_cube_size_, rngseeds_[levelmin_]);
}
store_rnd(levelmin_, randc[levelmin_]);
//... refinement levels
for (int ilevel = levelmin_ + 1; ilevel <= levelmax_; ++ilevel)
{
int lx[3], x0[3];
int shift[3], levelmin_poisson;
shift[0] = pcf_->getValue<int>("setup", "shift_x");
shift[1] = pcf_->getValue<int>("setup", "shift_y");
shift[2] = pcf_->getValue<int>("setup", "shift_z");
levelmin_poisson = pcf_->getValue<unsigned>("setup", "levelmin");
int lfac = 1 << (ilevel - levelmin_poisson);
std::array<int,3> margin;
if( prefh_->get_margin()>0 ){
margin[0] = prefh_->get_margin();
margin[1] = prefh_->get_margin();
margin[2] = prefh_->get_margin();
}else{
margin[0] = prefh_->size(ilevel, 0)/2;
margin[1] = prefh_->size(ilevel, 1)/2;
margin[2] = prefh_->size(ilevel, 2)/2;
}
lx[0] = prefh_->size(ilevel, 0) + 2 * margin[0];
lx[1] = prefh_->size(ilevel, 1) + 2 * margin[1];
lx[2] = prefh_->size(ilevel, 2) + 2 * margin[2];
x0[0] = prefh_->offset_abs(ilevel, 0) - lfac * shift[0] - margin[0];
x0[1] = prefh_->offset_abs(ilevel, 1) - lfac * shift[1] - margin[1];
x0[2] = prefh_->offset_abs(ilevel, 2) - lfac * shift[2] - margin[2];
if (randc[ilevel] == NULL)
randc[ilevel] =
new rng(*randc[ilevel - 1], ran_cube_size_, rngseeds_[ilevel], x0, lx);
delete randc[ilevel - 1];
randc[ilevel - 1] = NULL;
//... store numbers
store_rnd(ilevel, randc[ilevel]);
}
delete randc[levelmax_];
randc[levelmax_] = NULL;
}
void RNG_music::store_rnd(int ilevel, rng *prng)
{
int shift[3], levelmin_poisson;
shift[0] = pcf_->getValue<int>("setup", "shift_x");
shift[1] = pcf_->getValue<int>("setup", "shift_y");
shift[2] = pcf_->getValue<int>("setup", "shift_z");
levelmin_poisson = pcf_->getValue<unsigned>("setup", "levelmin");
int lfac = 1 << (ilevel - levelmin_poisson);
bool grafic_out = false;
if (grafic_out)
{
std::vector<float> data;
if (ilevel == levelmin_)
{
int N = 1 << levelmin_;
int i0, j0, k0;
i0 = -lfac * shift[0];
j0 = -lfac * shift[1];
k0 = -lfac * shift[2];
char fname[128];
sprintf(fname, "grafic_wnoise_%04d.bin", ilevel);
LOGUSER("Storing white noise field for grafic in file \'%s\'...", fname);
std::ofstream ofs(fname, std::ios::binary | std::ios::trunc);
data.assign(N * N, 0.0);
int blksize = 4 * sizeof(int);
int iseed = 0;
ofs.write(reinterpret_cast<char *>(&blksize), sizeof(int));
ofs.write(reinterpret_cast<char *>(&N), sizeof(int));
ofs.write(reinterpret_cast<char *>(&N), sizeof(int));
ofs.write(reinterpret_cast<char *>(&N), sizeof(int));
ofs.write(reinterpret_cast<char *>(&iseed), sizeof(int));
ofs.write(reinterpret_cast<char *>(&blksize), sizeof(int));
for (int k = 0; k < N; ++k)
{
#pragma omp parallel for
for (int j = 0; j < N; ++j)
for (int i = 0; i < N; ++i)
data[j * N + i] = -(*prng)(i + i0, j + j0, k + k0);
blksize = N * N * sizeof(float);
ofs.write(reinterpret_cast<char *>(&blksize), sizeof(int));
ofs.write(reinterpret_cast<char *>(&data[0]), N * N * sizeof(float));
ofs.write(reinterpret_cast<char *>(&blksize), sizeof(int));
}
ofs.close();
}
else
{
int nx, ny, nz;
int i0, j0, k0;
nx = prefh_->size(ilevel, 0);
ny = prefh_->size(ilevel, 1);
nz = prefh_->size(ilevel, 2);
i0 = prefh_->offset_abs(ilevel, 0) - lfac * shift[0];
j0 = prefh_->offset_abs(ilevel, 1) - lfac * shift[1];
k0 = prefh_->offset_abs(ilevel, 2) - lfac * shift[2];
char fname[128];
sprintf(fname, "grafic_wnoise_%04d.bin", ilevel);
LOGUSER("Storing white noise field for grafic in file \'%s\'...", fname);
LOGDEBUG("(%d,%d,%d) -- (%d,%d,%d) -- lfac = %d", nx, ny, nz, i0, j0, k0, lfac);
std::ofstream ofs(fname, std::ios::binary | std::ios::trunc);
data.assign(nx * ny, 0.0);
int blksize = 4 * sizeof(int);
int iseed = 0;
ofs.write(reinterpret_cast<char *>(&blksize), sizeof(int));
ofs.write(reinterpret_cast<char *>(&nz), sizeof(unsigned));
ofs.write(reinterpret_cast<char *>(&ny), sizeof(unsigned));
ofs.write(reinterpret_cast<char *>(&nx), sizeof(unsigned));
ofs.write(reinterpret_cast<char *>(&iseed), sizeof(int));
ofs.write(reinterpret_cast<char *>(&blksize), sizeof(int));
for (int k = 0; k < nz; ++k)
{
#pragma omp parallel for
for (int j = 0; j < ny; ++j)
for (int i = 0; i < nx; ++i)
data[j * nx + i] = -(*prng)(i + i0, j + j0, k + k0);
blksize = nx * ny * sizeof(float);
ofs.write(reinterpret_cast<char *>(&blksize), sizeof(int));
ofs.write(reinterpret_cast<char *>(&data[0]), nx * ny * sizeof(float));
ofs.write(reinterpret_cast<char *>(&blksize), sizeof(int));
}
ofs.close();
}
}
if (disk_cached_)
{
std::vector<real_t> data;
if (ilevel == levelmin_)
{
int N = 1 << levelmin_;
int i0, j0, k0;
i0 = -lfac * shift[0];
j0 = -lfac * shift[1];
k0 = -lfac * shift[2];
char fname[128];
sprintf(fname, "wnoise_%04d.bin", ilevel);
LOGUSER("Storing white noise field in file \'%s\'...", fname);
std::ofstream ofs(fname, std::ios::binary | std::ios::trunc);
ofs.write(reinterpret_cast<char *>(&N), sizeof(unsigned));
ofs.write(reinterpret_cast<char *>(&N), sizeof(unsigned));
ofs.write(reinterpret_cast<char *>(&N), sizeof(unsigned));
data.assign(N * N, 0.0);
for (int i = 0; i < N; ++i)
{
#pragma omp parallel for
for (int j = 0; j < N; ++j)
for (int k = 0; k < N; ++k)
data[j * N + k] = (*prng)(i + i0, j + j0, k + k0);
ofs.write(reinterpret_cast<char *>(&data[0]), N * N * sizeof(real_t));
}
ofs.close();
}
else
{
int nx, ny, nz;
int i0, j0, k0;
std::array<int,3> margin;
if( prefh_->get_margin()>0 ){
margin[0] = prefh_->get_margin();
margin[1] = prefh_->get_margin();
margin[2] = prefh_->get_margin();
}else{
margin[0] = prefh_->size(ilevel, 0)/2;
margin[1] = prefh_->size(ilevel, 1)/2;
margin[2] = prefh_->size(ilevel, 2)/2;
}
nx = prefh_->size(ilevel, 0) + 2 * margin[0];
ny = prefh_->size(ilevel, 1) + 2 * margin[1];
nz = prefh_->size(ilevel, 2) + 2 * margin[2];
i0 = prefh_->offset_abs(ilevel, 0) - lfac * shift[0] - margin[0];
j0 = prefh_->offset_abs(ilevel, 1) - lfac * shift[1] - margin[1];
k0 = prefh_->offset_abs(ilevel, 2) - lfac * shift[2] - margin[2];
char fname[128];
sprintf(fname, "wnoise_%04d.bin", ilevel);
LOGUSER("Storing white noise field in file \'%s\'...", fname);
std::ofstream ofs(fname, std::ios::binary | std::ios::trunc);
ofs.write(reinterpret_cast<char *>(&nx), sizeof(unsigned));
ofs.write(reinterpret_cast<char *>(&ny), sizeof(unsigned));
ofs.write(reinterpret_cast<char *>(&nz), sizeof(unsigned));
data.assign(ny * nz, 0.0);
for (int i = 0; i < nx; ++i)
{
#pragma omp parallel for
for (int j = 0; j < ny; ++j)
for (int k = 0; k < nz; ++k)
data[j * nz + k] = (*prng)(i + i0, j + j0, k + k0);
ofs.write(reinterpret_cast<char *>(&data[0]), ny * nz * sizeof(real_t));
}
ofs.close();
}
}
else
{
int nx, ny, nz;
int i0, j0, k0;
if (ilevel == levelmin_)
{
i0 = -lfac * shift[0];
j0 = -lfac * shift[1];
k0 = -lfac * shift[2];
nx = ny = nz = 1 << levelmin_;
}
else
{
std::array<int,3> margin;
if( prefh_->get_margin()>0 ){
margin[0] = prefh_->get_margin();
margin[1] = prefh_->get_margin();
margin[2] = prefh_->get_margin();
}else{
margin[0] = prefh_->size(ilevel, 0)/2;
margin[1] = prefh_->size(ilevel, 1)/2;
margin[2] = prefh_->size(ilevel, 2)/2;
}
nx = prefh_->size(ilevel, 0) + 2 * margin[0];
ny = prefh_->size(ilevel, 1) + 2 * margin[1];
nz = prefh_->size(ilevel, 2) + 2 * margin[2];
i0 = prefh_->offset_abs(ilevel, 0) - lfac * shift[0] - margin[0];
j0 = prefh_->offset_abs(ilevel, 1) - lfac * shift[1] - margin[1];
k0 = prefh_->offset_abs(ilevel, 2) - lfac * shift[2] - margin[2];
}
mem_cache_[ilevel - levelmin_] = new std::vector<real_t>(nx * ny * nz, 0.0);
LOGUSER("Copying white noise to mem cache...");
#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)
(*mem_cache_[ilevel - levelmin_])[((size_t)i * ny + (size_t)j) * nz + (size_t)k] =
(*prng)(i + i0, j + j0, k + k0);
}
}
void RNG_music::fill_grid(int ilevel, DensityGrid<real_t> &A)
{
if (!initialized_)
{
LOGERR("Call to RNG_music::fill_grid before call to RNG_music::initialize_for_grid_structure");
throw std::runtime_error("invalid call order for random number generator");
}
if (restart_)
LOGINFO("Attempting to restart using random numbers for level %d\n from file \'wnoise_%04d.bin\'.", ilevel,
ilevel);
if (disk_cached_)
{
char fname[128];
sprintf(fname, "wnoise_%04d.bin", ilevel);
LOGUSER("Loading white noise from file \'%s\'...", fname);
std::ifstream ifs(fname, std::ios::binary);
if (!ifs.good())
{
LOGERR("White noise file \'%s\'was not found.", fname);
throw std::runtime_error("A white noise file was not found. This is an internal inconsistency and bad.");
}
int nx, ny, nz;
ifs.read(reinterpret_cast<char *>(&nx), sizeof(int));
ifs.read(reinterpret_cast<char *>(&ny), sizeof(int));
ifs.read(reinterpret_cast<char *>(&nz), sizeof(int));
if (nx != (int)A.size(0) || ny != (int)A.size(1) || nz != (int)A.size(2))
{
std::array<int,3> margin;
if( prefh_->get_margin()>0 ){
margin[0] = prefh_->get_margin();
margin[1] = prefh_->get_margin();
margin[2] = prefh_->get_margin();
}else{
margin[0] = prefh_->size(ilevel, 0)/2;
margin[1] = prefh_->size(ilevel, 1)/2;
margin[2] = prefh_->size(ilevel, 2)/2;
}
if (nx == (int)A.size(0) + 2 * margin[0] && ny == (int)A.size(1) + 2 * margin[1] && nz == (int)A.size(2) + 2 * margin[2])
{
int ox = margin[0], oy = margin[1], oz = margin[2];
std::vector<real_t> slice(ny * nz, 0.0);
for (int i = 0; i < nx; ++i)
{
ifs.read(reinterpret_cast<char *>(&slice[0]), ny * nz * sizeof(real_t));
if (i < ox)
continue;
if (i >= 3 * ox)
break;
#pragma omp parallel for
for (int j = oy; j < 3 * oy; ++j)
for (int k = oz; k < 3 * oz; ++k)
A(i - ox, j - oy, k - oz) = slice[j * nz + k];
}
ifs.close();
}
else
{
LOGERR("White noise file is not aligned with array. File: [%d,%d,%d]. Mem: [%d,%d,%d].", nx, ny, nz, A.size(0),
A.size(1), A.size(2));
throw std::runtime_error(
"White noise file is not aligned with array. This is an internal inconsistency and bad.");
}
}
else
{
for (int i = 0; i < nx; ++i)
{
std::vector<real_t> slice(ny * nz, 0.0);
ifs.read(reinterpret_cast<char *>(&slice[0]), ny * nz * sizeof(real_t));
#pragma omp parallel for
for (int j = 0; j < ny; ++j)
for (int k = 0; k < nz; ++k)
A(i, j, k) = slice[j * nz + k];
}
ifs.close();
}
}
else
{
LOGUSER("Copying white noise from memory cache...");
if (mem_cache_[ilevel - levelmin_] == NULL)
LOGERR("Tried to access mem-cached random numbers for level %d. But these are not available!\n", ilevel);
int nx(A.size(0)), ny(A.size(1)), nz(A.size(2));
if ((size_t)nx * (size_t)ny * (size_t)nz != mem_cache_[ilevel - levelmin_]->size())
{
LOGERR("White noise file is not aligned with array. File: [%d,%d,%d]. Mem: [%d,%d,%d].", nx, ny, nz, A.size(0),
A.size(1), A.size(2));
throw std::runtime_error("White noise file is not aligned with array. This is an internal inconsistency and bad");
}
#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)
A(i, j, k) = (*mem_cache_[ilevel - levelmin_])[((size_t)i * ny + (size_t)j) * nz + (size_t)k];
std::vector<real_t>().swap(*mem_cache_[ilevel - levelmin_]);
delete mem_cache_[ilevel - levelmin_];
mem_cache_[ilevel - levelmin_] = NULL;
}
};
namespace
{
RNG_plugin_creator_concrete<RNG_music> creator("MUSIC");
}

1030
src/plugins/random_music_wnoise_generator.cc Normal file
View file

File diff suppressed because it is too large Load diff

197
src/plugins/random_music_wnoise_generator.hh Normal file
View file

@ -0,0 +1,197 @@
#pragma once
#include <vector>
#include <map>
#include "general.hh"
#include "mesh.hh"
#define DEF_RAN_CUBE_SIZE 32
/*!
* @brief encapsulates all things random number generator related
*/
template <typename T>
class music_wnoise_generator
{
public:
unsigned
res_, //!< resolution of the full mesh
cubesize_, //!< size of one independent random number cube
ncubes_; //!< number of random number cubes to cover the full mesh
long baseseed_; //!< base seed from which cube seeds are computed
protected:
//! vector of 3D meshes (the random number cubes) with random numbers
std::vector<Meshvar<T> *> rnums_;
//! map of 3D indices to cube index
std::map<size_t, size_t> cubemap_;
typedef std::map<size_t, size_t>::iterator cubemap_iterator;
protected:
//! register a cube with the hash map
void register_cube(int i, int j, int k);
//! fills a subcube with random numbers
double fill_cube(int i, int j, int k);
//! subtract a constant from an entire cube
void subtract_from_cube(int i, int j, int k, double val);
//! copy random numbers from a cube to a full grid array
template <class C>
void copy_cube(int i, int j, int k, C &dat)
{
int offi, offj, offk;
offi = i * cubesize_;
offj = j * cubesize_;
offk = k * cubesize_;
i = (i + ncubes_) % ncubes_;
j = (j + ncubes_) % ncubes_;
k = (k + ncubes_) % ncubes_;
size_t icube = (i * ncubes_ + j) * ncubes_ + k;
cubemap_iterator it = cubemap_.find(icube);
if (it == cubemap_.end())
{
LOGERR("attempting to copy data from non-existing RND cube %d,%d,%d", i, j, k);
throw std::runtime_error("attempting to copy data from non-existing RND cube");
}
size_t cubeidx = it->second;
for (int ii = 0; ii < (int)cubesize_; ++ii)
for (int jj = 0; jj < (int)cubesize_; ++jj)
for (int kk = 0; kk < (int)cubesize_; ++kk)
dat(offi + ii, offj + jj, offk + kk) = (*rnums_[cubeidx])(ii, jj, kk);
}
//! free the memory associated with a subcube
void free_cube(int i, int j, int k);
//! initialize member variables and allocate memory
void initialize(void);
//! fill a cubic subvolume of the full grid with random numbers
double fill_subvolume(int *i0, int *n);
//! fill an entire grid with random numbers
double fill_all(void);
//! fill an external array instead of the internal field
template <class C>
double fill_all(C &dat)
{
double sum = 0.0;
for (int i = 0; i < (int)ncubes_; ++i)
for (int j = 0; j < (int)ncubes_; ++j)
for (int k = 0; k < (int)ncubes_; ++k)
{
int ii(i), jj(j), kk(k);
register_cube(ii, jj, kk);
}
#pragma omp parallel for reduction(+ \
: sum)
for (int i = 0; i < (int)ncubes_; ++i)
for (int j = 0; j < (int)ncubes_; ++j)
for (int k = 0; k < (int)ncubes_; ++k)
{
int ii(i), jj(j), kk(k);
ii = (ii + ncubes_) % ncubes_;
jj = (jj + ncubes_) % ncubes_;
kk = (kk + ncubes_) % ncubes_;
sum += fill_cube(ii, jj, kk);
copy_cube(ii, jj, kk, dat);
free_cube(ii, jj, kk);
}
return sum / (ncubes_ * ncubes_ * ncubes_);
}
//! write the number of allocated random number cubes to stdout
void print_allocated(void);
public:
//! constructor
music_wnoise_generator(unsigned res, unsigned cubesize, long baseseed, int *x0, int *lx);
//! constructor for constrained fine field
music_wnoise_generator(music_wnoise_generator<T> &rc, unsigned cubesize, long baseseed, int *x0_ = NULL, int *lx_ = NULL, bool zeromean = true);
//! constructor
music_wnoise_generator(unsigned res, unsigned cubesize, long baseseed, bool zeromean = true);
//! constructor to read white noise from file
music_wnoise_generator(unsigned res, std::string randfname, bool rndsign);
//! copy constructor for averaged field (not copying) hence explicit!
explicit music_wnoise_generator(/*const*/ music_wnoise_generator<T> &rc);
//! destructor
~music_wnoise_generator()
{
for (unsigned i = 0; i < rnums_.size(); ++i)
if (rnums_[i] != NULL)
delete rnums_[i];
rnums_.clear();
}
//! access a random number, this allocates a cube and fills it with consistent random numbers
inline T &operator()(int i, int j, int k, bool fillrand = true)
{
int ic, jc, kc, is, js, ks;
if (ncubes_ == 0)
throw std::runtime_error("music_wnoise_generator: internal error, not properly initialized");
//... determine cube
ic = (int)((double)i / cubesize_ + ncubes_) % ncubes_;
jc = (int)((double)j / cubesize_ + ncubes_) % ncubes_;
kc = (int)((double)k / cubesize_ + ncubes_) % ncubes_;
size_t icube = ((size_t)ic * ncubes_ + (size_t)jc) * ncubes_ + (size_t)kc;
cubemap_iterator it = cubemap_.find(icube);
if (it == cubemap_.end())
{
LOGERR("Attempting to copy data from non-existing RND cube %d,%d,%d @ %d,%d,%d", ic, jc, kc, i, j, k);
throw std::runtime_error("attempting to copy data from non-existing RND cube");
}
size_t cubeidx = it->second;
if (rnums_[cubeidx] == NULL)
{
LOGERR("Attempting to access data from non-allocated RND cube %d,%d,%d", ic, jc, kc);
throw std::runtime_error("attempting to access data from non-allocated RND cube");
}
//... determine cell in cube
is = (i - ic * cubesize_ + cubesize_) % cubesize_;
js = (j - jc * cubesize_ + cubesize_) % cubesize_;
ks = (k - kc * cubesize_ + cubesize_) % cubesize_;
return (*rnums_[cubeidx])(is, js, ks);
}
//! free all cubes
void free_all_mem(void)
{
for (unsigned i = 0; i < rnums_.size(); ++i)
if (rnums_[i] != NULL)
{
delete rnums_[i];
rnums_[i] = NULL;
}
}
};

887
src/plugins/random_panphasia.cc Normal file
View file

@ -0,0 +1,887 @@
#ifdef HAVE_PANPHASIA
#include "random.hh"
#include <cctype>
#include <cstring>
#include <stdint.h>
#include "densities.hh"
#include "HDF_IO.hh"
const int maxdim = 60, maxlev = 50, maxpow = 3 * maxdim;
typedef int rand_offset_[5];
typedef struct
{
int state[133]; // Nstore = Nstate (=5) + Nbatch (=128)
int need_fill;
int pos;
} rand_state_;
/* pan_state_ struct -- corresponds to respective fortran module in panphasia_routines.f
* data structure that contains all panphasia state variables
* it needs to get passed between the fortran routines to enable
* thread-safe execution.
*/
typedef struct
{
int base_state[5], base_lev_start[5][maxdim + 1];
rand_offset_ poweroffset[maxpow + 1], superjump;
rand_state_ current_state[maxpow + 2];
int layer_min, layer_max, indep_field;
long long xorigin_store[2][2][2], yorigin_store[2][2][2], zorigin_store[2][2][2];
int lev_common, layer_min_store, layer_max_store;
long long ix_abs_store, iy_abs_store, iz_abs_store, ix_per_store, iy_per_store, iz_per_store, ix_rel_store,
iy_rel_store, iz_rel_store;
double exp_coeffs[8][8][maxdim + 2];
long long xcursor[maxdim + 1], ycursor[maxdim + 1], zcursor[maxdim + 1];
int ixshift[2][2][2], iyshift[2][2][2], izshift[2][2][2];
double cell_data[9][8];
int ixh_last, iyh_last, izh_last;
int init;
int init_cell_props;
int init_lecuyer_state;
long long p_xcursor[62], p_ycursor[62], p_zcursor[62];
} pan_state_;
extern "C"
{
void start_panphasia_(pan_state_ *lstate, const char *descriptor, int *ngrid, int *bverbose);
void parse_descriptor_(const char *descriptor, int16_t *l, int32_t *ix, int32_t *iy, int32_t *iz, int16_t *side1,
int16_t *side2, int16_t *side3, int32_t *check_int, char *name);
void panphasia_cell_properties_(pan_state_ *lstate, int *ixcell, int *iycell, int *izcell, double *cell_prop);
void adv_panphasia_cell_properties_(pan_state_ *lstate, int *ixcell, int *iycell, int *izcell, int *layer_min,
int *layer_max, int *indep_field, double *cell_prop);
void set_phases_and_rel_origin_(pan_state_ *lstate, const char *descriptor, int *lev, long long *ix_rel,
long long *iy_rel, long long *iz_rel, int *VERBOSE);
/*void set_local_box_( pan_state_ *lstate, int lev, int8_t ix_abs, int8_t iy_abs, int8_t iz_abs,
int8_t ix_per, int8_t iy_per, int8_t iz_per, int8_t ix_rel, int8_t iy_rel,
int8_t iz_rel, int wn_level_base, int8_t check_rand, char *phase_name, int MYID);*/
/*extern struct {
int layer_min, layer_max, hoswitch;
}oct_range_;
*/
}
class RNG_panphasia : public RNG_plugin
{
private:
void forward_transform_field(real_t *field, int n0, int n1, int n2);
void forward_transform_field(real_t *field, int n) { forward_transform_field(field, n, n, n); }
void backward_transform_field(real_t *field, int n0, int n1, int n2);
void backward_transform_field(real_t *field, int n) { backward_transform_field(field, n, n, n); }
protected:
std::string descriptor_string_;
int num_threads_;
int levelmin_, levelmin_final_, levelmax_, ngrid_;
bool incongruent_fields_;
double inter_grid_phase_adjustment_;
// double translation_phase_;
pan_state_ *lstate;
int grid_p_, grid_m_;
double grid_rescale_fac_;
int coordinate_system_shift_[3];
int ix_abs_[3], ix_per_[3], ix_rel_[3], level_p_, lextra_;
const refinement_hierarchy *prefh_;
std::array<int,3> margins_;
struct panphasia_descriptor
{
int16_t wn_level_base;
int32_t i_xorigin_base, i_yorigin_base, i_zorigin_base;
int16_t i_base, i_base_y, i_base_z;
int32_t check_rand;
std::string name;
explicit panphasia_descriptor(std::string dstring)
{
char tmp[100];
memset(tmp, ' ', 100);
parse_descriptor_(dstring.c_str(), &wn_level_base, &i_xorigin_base, &i_yorigin_base, &i_zorigin_base, &i_base,
&i_base_y, &i_base_z, &check_rand, tmp);
for (int i = 0; i < 100; i++)
if (tmp[i] == ' ')
{
tmp[i] = '\0';
break;
}
name = tmp;
name.erase(std::remove(name.begin(), name.end(), ' '), name.end());
}
};
void clear_panphasia_thread_states(void)
{
for (int i = 0; i < num_threads_; ++i)
{
lstate[i].init = 0;
lstate[i].init_cell_props = 0;
lstate[i].init_lecuyer_state = 0;
}
}
// greatest common divisor
int gcd(int a, int b)
{
if (b == 0)
return a;
return gcd(b, a % b);
}
// least common multiple
int lcm(int a, int b) { return abs(a * b) / gcd(a, b); }
// Two or largest power of 2 less than the argument
int largest_power_two_lte(int b)
{
int a = 1;
if (b <= a)
return a;
while (2 * a < b)
a = 2 * a;
return a;
}
panphasia_descriptor *pdescriptor_;
public:
explicit RNG_panphasia(config_file &cf) : RNG_plugin(cf)
{
descriptor_string_ = pcf_->getValue<std::string>("random", "descriptor");
#ifdef _OPENMP
num_threads_ = omp_get_max_threads();
#else
num_threads_ = 1;
#endif
// create independent state descriptions for each thread
lstate = new pan_state_[num_threads_];
// parse the descriptor for its properties
pdescriptor_ = new panphasia_descriptor(descriptor_string_);
LOGINFO("PANPHASIA: descriptor \'%s\' is base %d,", pdescriptor_->name.c_str(), pdescriptor_->i_base);
// write panphasia base size into config file for the grid construction
// as the gridding unit we use the least common multiple of 2 and i_base
std::stringstream ss;
// ARJ ss << lcm(2, pdescriptor_->i_base);
// ss << two_or_largest_power_two_less_than(pdescriptor_->i_base);//ARJ
ss << 2; // ARJ - set gridding unit to two
pcf_->insertValue("setup", "gridding_unit", ss.str());
ss.str(std::string());
ss << pdescriptor_->i_base;
pcf_->insertValue("random", "base_unit", ss.str());
pcf_->insertValue("setup","fourier_splicing","false");
}
void initialize_for_grid_structure(const refinement_hierarchy &refh)
{
prefh_ = &refh;
levelmin_ = prefh_->levelmin();
levelmin_final_ = pcf_->getValue<unsigned>("setup", "levelmin");
levelmax_ = prefh_->levelmax();
if( refh.get_margin() < 0 ){
margins_ = {-1,-1,-1};
}else{
margins_ = { refh.get_margin(), refh.get_margin(), refh.get_margin() };
}
clear_panphasia_thread_states();
LOGINFO("PANPHASIA: running with %d threads", num_threads_);
// if ngrid is not a multiple of i_base, then we need to enlarge and then sample down
ngrid_ = 1 << levelmin_;
grid_p_ = pdescriptor_->i_base;
grid_m_ = largest_power_two_lte(grid_p_);
lextra_ = (log10((double)ngrid_ / (double)pdescriptor_->i_base) + 0.001) / log10(2.0);
int ratio = 1 << lextra_;
grid_rescale_fac_ = 1.0;
coordinate_system_shift_[0] = -pcf_->getValue<int>("setup", "shift_x");
coordinate_system_shift_[1] = -pcf_->getValue<int>("setup", "shift_y");
coordinate_system_shift_[2] = -pcf_->getValue<int>("setup", "shift_z");
incongruent_fields_ = false;
if (ngrid_ != ratio * pdescriptor_->i_base)
{
incongruent_fields_ = true;
ngrid_ = 2 * ratio * pdescriptor_->i_base;
grid_rescale_fac_ = (double)ngrid_ / (1 << levelmin_);
LOGINFO("PANPHASIA: will use a higher resolution:\n"
" (%d -> %d) * 2**ref compatible with PANPHASIA\n"
" will Fourier interpolate after.",
grid_m_, grid_p_);
}
}
~RNG_panphasia() { delete[] lstate; }
void fill_grid(int level, DensityGrid<real_t> &R);
bool is_multiscale() const { return true; }
};
void RNG_panphasia::forward_transform_field(real_t *field, int nx, int ny, int nz)
{
fftw_real *rfield = reinterpret_cast<fftw_real *>(field);
fftw_complex *cfield = reinterpret_cast<fftw_complex *>(field);
#ifdef FFTW3
#ifdef SINGLE_PRECISION
fftwf_plan pf = fftwf_plan_dft_r2c_3d(nx, ny, nz, rfield, cfield, FFTW_ESTIMATE);
#else
fftw_plan pf = fftw_plan_dft_r2c_3d(nx, ny, nz, rfield, cfield, FFTW_ESTIMATE);
#endif
#else
rfftwnd_plan pf = rfftw3d_create_plan(nx, ny, nz, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE | FFTW_IN_PLACE);
#endif
#ifdef FFTW3
#ifdef SINGLE_PRECISION
fftwf_execute(pf);
#else
fftw_execute(pf);
#endif
#else
#ifndef SINGLETHREAD_FFTW
rfftwnd_threads_one_real_to_complex(num_threads_, pf, rfield, NULL);
#else
rfftwnd_one_real_to_complex(pf, rfield, NULL);
#endif
#endif
}
void RNG_panphasia::backward_transform_field(real_t *field, int nx, int ny, int nz)
{
fftw_real *rfield = reinterpret_cast<fftw_real *>(field);
fftw_complex *cfield = reinterpret_cast<fftw_complex *>(field);
#ifdef FFTW3
#ifdef SINGLE_PRECISION
fftwf_plan ipf = fftwf_plan_dft_c2r_3d(nx, ny, nz, cfield, rfield, FFTW_ESTIMATE);
#else
fftw_plan ipf = fftw_plan_dft_c2r_3d(nx, ny, nz, cfield, rfield, FFTW_ESTIMATE);
#endif
#else
rfftwnd_plan ipf = rfftw3d_create_plan(nx, ny, nz, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE | FFTW_IN_PLACE);
#endif
#ifdef FFTW3
#ifdef SINGLE_PRECISION
fftwf_execute(ipf);
#else
fftw_execute(ipf);
#endif
#else
#ifndef SINGLETHREAD_FFTW
rfftwnd_threads_one_complex_to_real(num_threads_, ipf, cfield, NULL);
#else
rfftwnd_one_complex_to_real(ipf, cfield, NULL);
#endif
#endif
}
#include <sys/time.h>
inline double get_wtime(void)
{
#ifdef _OPENMP
return omp_get_wtime();
#else
return (double)clock() / CLOCKS_PER_SEC;
#endif
}
void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
{
fftw_real *pr0, *pr1, *pr2, *pr3, *pr4;
fftw_complex *pc0, *pc1, *pc2, *pc3, *pc4;
// determine resolution and offset so that we can do proper resampling
int ileft[3], ileft_corner[3], nx[3], nxremap[3];
int iexpand_left[3];
for (int k = 0; k < 3; ++k)
{
ileft[k] = prefh_->offset_abs(level, k);
nx[k] = R.size(k);
assert(nx[k] % 4 == 0);
if (level == levelmin_)
{
ileft_corner[k] = ileft[k]; // Top level - periodic
}
else
{
if( margins_[0] < 0 ){
ileft_corner[k] = (ileft[k] - nx[k] / 4 + (1 << level)) % (1 << level); // Isolated
ileft_corner[k] = (ileft[k] - nx[k] / 4 + (1 << level)) % (1 << level); // Isolated
}else{
ileft_corner[k] = (ileft[k] - margins_[k] + (1 << level)) % (1 << level); // Isolated
ileft_corner[k] = (ileft[k] - margins_[k] + (1 << level)) % (1 << level); // Isolated
}
}
iexpand_left[k] = (ileft_corner[k] % grid_m_ == 0) ? 0 : ileft_corner[k] % grid_m_;
// fprintf(stderr, "dim=%c : ileft = %d, ileft_corner %d, nx = %d\n", 'x' + k, ileft[k],ileft_corner[k],nx[k]);
};
int ileft_corner_m[3], ileft_corner_p[3], nx_m[3];
int ileft_max_expand = std::max(iexpand_left[0], std::max(iexpand_left[1], iexpand_left[2]));
for (int k = 0; k < 3; ++k)
{
ileft_corner_m[k] = ((ileft_corner[k] - iexpand_left[k]) +
coordinate_system_shift_[k] * (1 << (level - levelmin_final_)) + (1 << level)) %
(1 << level);
ileft_corner_p[k] = grid_p_ * ileft_corner_m[k] / grid_m_;
nx_m[k] = (ileft_max_expand != 0) ? nx[k] + ileft_max_expand : nx[k];
if (nx_m[k] % grid_m_ != 0)
nx_m[k] = nx_m[k] + grid_m_ - nx_m[k] % grid_m_;
nxremap[k] = grid_p_ * nx_m[k] / grid_m_;
if (nxremap[k] % 2 == 1)
{
nx_m[k] = nx_m[k] + grid_m_;
nxremap[k] = grid_p_ * nx_m[k] / grid_m_;
}
}
if ((nx_m[0] != nx_m[1]) || (nx_m[0] != nx_m[2]))
LOGERR("Fatal error: non-cubic refinement being requested");
inter_grid_phase_adjustment_ = M_PI * (1.0 / (double)nx_m[0] - 1.0 / (double)nxremap[0]);
// LOGINFO("The value of the phase adjustement is %f\n", inter_grid_phase_adjustment_);
// LOGINFO("ileft[0],ileft[1],ileft[2] %d %d %d", ileft[0], ileft[1], ileft[2]);
// LOGINFO("ileft_corner[0,1,2] %d %d %d", ileft_corner[0], ileft_corner[1], ileft_corner[2]);
// LOGINFO("iexpand_left[1,2,3] = (%d, %d, %d) Max %d ",iexpand_left[0],iexpand_left[1],iexpand_left[2], ileft_max_expand);
// LOGINFO("ileft_corner_m[0,1,2] = (%d,%d,%d)",ileft_corner_m[0],ileft_corner_m[1],ileft_corner_m[2]);
// LOGINFO("grid_m_ %d grid_p_ %d",grid_m_,grid_p_);
// LOGINFO("nx_m[0,1,2] = (%d,%d,%d)",nx_m[0],nx_m[1],nx_m[2]);
// LOGINFO("ileft_corner_p[0,1,2] = (%d,%d,%d)",ileft_corner_p[0],ileft_corner_p[1],ileft_corner_p[2]);
// LOGINFO("nxremap[0,1,2] = (%d,%d,%d)",nxremap[0],nxremap[1],nxremap[2]);
size_t ngp = nxremap[0] * nxremap[1] * (nxremap[2] + 2);
pr0 = new fftw_real[ngp];
pr1 = new fftw_real[ngp];
pr2 = new fftw_real[ngp];
pr3 = new fftw_real[ngp];
pr4 = new fftw_real[ngp];
pc0 = reinterpret_cast<fftw_complex *>(pr0);
pc1 = reinterpret_cast<fftw_complex *>(pr1);
pc2 = reinterpret_cast<fftw_complex *>(pr2);
pc3 = reinterpret_cast<fftw_complex *>(pr3);
pc4 = reinterpret_cast<fftw_complex *>(pr4);
LOGINFO("calculating PANPHASIA random numbers for level %d...", level);
clear_panphasia_thread_states();
double t1 = get_wtime();
double tp = t1;
#pragma omp parallel
{
#ifdef _OPENMP
const int mythread = omp_get_thread_num();
#else
const int mythread = 0;
#endif
int odd_x, odd_y, odd_z;
int ng_level = ngrid_ * (1 << (level - levelmin_)); // full resolution of current level
int verbosity = (mythread == 0);
char descriptor[100];
memset(descriptor, 0, 100);
memcpy(descriptor, descriptor_string_.c_str(), descriptor_string_.size());
if (level == levelmin_)
{
start_panphasia_(&lstate[mythread], descriptor, &ng_level, &verbosity);
}
{
int level_p, lextra;
long long ix_rel[3];
panphasia_descriptor d(descriptor_string_);
lextra = (log10((double)ng_level / (double)d.i_base) + 0.001) / log10(2.0);
level_p = d.wn_level_base + lextra;
int ratio = 1 << lextra;
assert(ng_level == ratio * d.i_base);
ix_rel[0] = ileft_corner_p[0];
ix_rel[1] = ileft_corner_p[1];
ix_rel[2] = ileft_corner_p[2];
// Code above ignores the coordinate_system_shift_ - but currently this is set to zero //
lstate[mythread].layer_min = 0;
lstate[mythread].layer_max = level_p;
lstate[mythread].indep_field = 1;
set_phases_and_rel_origin_(&lstate[mythread], descriptor, &level_p, &ix_rel[0], &ix_rel[1], &ix_rel[2],
&verbosity);
// LOGUSER(" called set_phases_and_rel_origin level %d ix_rel iy_rel iz_rel %d %d %d\n", level_p, ix_rel[0], ix_rel[1], ix_rel[2]);
odd_x = ix_rel[0] % 2;
odd_y = ix_rel[1] % 2;
odd_z = ix_rel[2] % 2;
}
if (verbosity)
t1 = get_wtime();
//***************************************************************
// Process Panphasia values: p000, p001, p010, p100 and indep field
//****************************************************************
// START //
#pragma omp for // nowait
for (int i = 0; i < nxremap[0] / 2 + odd_x; ++i)
{
double cell_prop[9];
pan_state_ *ps = &lstate[mythread];
for (int j = 0; j < nxremap[1] / 2 + odd_y; ++j)
for (int k = 0; k < nxremap[2] / 2 + odd_z; ++k)
{
// ARJ - added inner set of loops to speed up evaluation of Panphasia
for (int ix = 0; ix < 2; ++ix)
for (int iy = 0; iy < 2; ++iy)
for (int iz = 0; iz < 2; ++iz)
{
int ii = 2 * i + ix - odd_x;
int jj = 2 * j + iy - odd_y;
int kk = 2 * k + iz - odd_z;
if (((ii >= 0) && (ii < nxremap[0])) && ((jj >= 0) && (jj < nxremap[1])) &&
((kk >= 0) && (kk < nxremap[2])))
{
size_t idx = ((size_t)ii * nxremap[1] + (size_t)jj) * (nxremap[2] + 2) + (size_t)kk;
adv_panphasia_cell_properties_(ps, &ii, &jj, &kk, &ps->layer_min, &ps->layer_max, &ps->indep_field,
cell_prop);
pr0[idx] = cell_prop[0];
pr1[idx] = cell_prop[4];
pr2[idx] = cell_prop[2];
pr3[idx] = cell_prop[1];
pr4[idx] = cell_prop[8];
}
}
}
}
}
LOGUSER("time for calculating PANPHASIA for level %d : %f s, %f µs/cell", level, get_wtime() - t1,
1e6 * (get_wtime() - t1) / ((double)nxremap[2] * (double)nxremap[1] * (double)nxremap[0]));
LOGUSER("time for calculating PANPHASIA for level %d : %f s, %f µs/cell", level, get_wtime() - t1,
1e6 * (get_wtime() - t1) / ((double)nxremap[2] * (double)nxremap[1] * (double)nxremap[0]));
//////////////////////////////////////////////////////////////////////////////////////////////
LOGUSER("\033[31mtiming level %d [adv_panphasia_cell_properties]: %f s\033[0m", level, get_wtime() - tp);
tp = get_wtime();
/////////////////////////////////////////////////////////////////////////
// transform and convolve with Legendres
forward_transform_field(pr0, nxremap[0], nxremap[1], nxremap[2]);
forward_transform_field(pr1, nxremap[0], nxremap[1], nxremap[2]);
forward_transform_field(pr2, nxremap[0], nxremap[1], nxremap[2]);
forward_transform_field(pr3, nxremap[0], nxremap[1], nxremap[2]);
forward_transform_field(pr4, nxremap[0], nxremap[1], nxremap[2]);
#pragma omp parallel for
for (int i = 0; i < nxremap[0]; i++)
for (int j = 0; j < nxremap[1]; j++)
for (int k = 0; k < nxremap[2] / 2 + 1; k++)
{
size_t idx = ((size_t)i * nxremap[1] + (size_t)j) * (nxremap[2] / 2 + 1) + (size_t)k;
double fx(1.0), fy(1.0), fz(1.0), arg = 0.;
complex gx(0., 0.), gy(0., 0.), gz(0., 0.);
int ii(i), jj(j), kk(k);
if (i > nxremap[0] / 2)
ii -= nxremap[0];
if (j > nxremap[1] / 2)
jj -= nxremap[1];
// int kkmax = std::max(abs(ii),std::max(abs(jj),abs(kk)));
if (ii != 0)
{
arg = M_PI * (double)ii / (double)nxremap[0];
fx = sin(arg) / arg;
gx = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
}
else
{
fx = 1.0;
gx = 0.0;
}
if (jj != 0)
{
arg = M_PI * (double)jj / (double)nxremap[1];
fy = sin(arg) / arg;
gy = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
}
else
{
fy = 1.0;
gy = 0.0;
}
if (kk != 0)
{
arg = M_PI * (double)kk / (double)nxremap[2];
fz = sin(arg) / arg;
gz = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
}
else
{
fz = 1.0;
gz = 0.0;
}
complex temp_comp = (fx + sqrt(3.0) * gx) * (fy + sqrt(3.0) * gy) * (fz + sqrt(3.0) * gz);
double magnitude = sqrt(1.0 - std::abs(temp_comp * temp_comp));
if (abs(ii) != nxremap[0] / 2 && abs(jj) != nxremap[1] / 2 &&
abs(kk) != nxremap[2] / 2)
{ // kkmax != nxremap[2]/2 ){
complex x, y0(RE(pc0[idx]), IM(pc0[idx])), y1(RE(pc1[idx]), IM(pc1[idx])), y2(RE(pc2[idx]), IM(pc2[idx])),
y3(RE(pc3[idx]), IM(pc3[idx])), y4(RE(pc4[idx]), IM(pc4[idx]));
x = y0 * fx * fy * fz + sqrt(3.0) * (y1 * gx * fy * fz + y2 * fx * gy * fz + y3 * fx * fy * gz) +
y4 * magnitude;
RE(pc0[idx]) = x.real();
IM(pc0[idx]) = x.imag();
}
}
// END
LOGUSER("\033[31mtiming level %d [build panphasia field]: %f s\033[0m", level, get_wtime() - tp);
tp = get_wtime();
//***************************************************************
// Process Panphasia values: p000, p001, p010, p100 and indep field
//****************************************************************
#pragma omp parallel
{
#ifdef _OPENMP
const int mythread = omp_get_thread_num();
#else
const int mythread = 0;
#endif
int odd_x, odd_y, odd_z;
int ng_level = ngrid_ * (1 << (level - levelmin_)); // full resolution of current level
int verbosity = (mythread == 0);
char descriptor[100];
memset(descriptor, 0, 100);
memcpy(descriptor, descriptor_string_.c_str(), descriptor_string_.size());
if (level == levelmin_)
{
start_panphasia_(&lstate[mythread], descriptor, &ng_level, &verbosity);
}
{
int level_p, lextra;
long long ix_rel[3];
panphasia_descriptor d(descriptor_string_);
lextra = (log10((double)ng_level / (double)d.i_base) + 0.001) / log10(2.0);
level_p = d.wn_level_base + lextra;
int ratio = 1 << lextra;
assert(ng_level == ratio * d.i_base);
ix_rel[0] = ileft_corner_p[0];
ix_rel[1] = ileft_corner_p[1];
ix_rel[2] = ileft_corner_p[2];
// Code above ignores the coordinate_system_shift_ - but currently this is set to zero //
lstate[mythread].layer_min = 0;
lstate[mythread].layer_max = level_p;
lstate[mythread].indep_field = 1;
set_phases_and_rel_origin_(&lstate[mythread], descriptor, &level_p, &ix_rel[0], &ix_rel[1], &ix_rel[2],
&verbosity);
// LOGINFO(" called set_phases_and_rel_origin level %d ix_rel iy_rel iz_rel %d %d %d\n", level_p, ix_rel[0], ix_rel[1], ix_rel[2]);
odd_x = ix_rel[0] % 2;
odd_y = ix_rel[1] % 2;
odd_z = ix_rel[2] % 2;
}
if (verbosity)
t1 = get_wtime();
// START //
//***************************************************************
// Process Panphasia values: p110, p011, p101, p111
//****************************************************************
#pragma omp for // nowait
for (int i = 0; i < nxremap[0] / 2 + odd_x; ++i)
{
double cell_prop[9];
pan_state_ *ps = &lstate[mythread];
for (int j = 0; j < nxremap[1] / 2 + odd_y; ++j)
for (int k = 0; k < nxremap[2] / 2 + odd_z; ++k)
{
// ARJ - added inner set of loops to speed up evaluation of Panphasia
for (int ix = 0; ix < 2; ++ix)
for (int iy = 0; iy < 2; ++iy)
for (int iz = 0; iz < 2; ++iz)
{
int ii = 2 * i + ix - odd_x;
int jj = 2 * j + iy - odd_y;
int kk = 2 * k + iz - odd_z;
if (((ii >= 0) && (ii < nxremap[0])) && ((jj >= 0) && (jj < nxremap[1])) &&
((kk >= 0) && (kk < nxremap[2])))
{
size_t idx = ((size_t)ii * nxremap[1] + (size_t)jj) * (nxremap[2] + 2) + (size_t)kk;
adv_panphasia_cell_properties_(ps, &ii, &jj, &kk, &ps->layer_min, &ps->layer_max, &ps->indep_field,
cell_prop);
pr1[idx] = cell_prop[6];
pr2[idx] = cell_prop[3];
pr3[idx] = cell_prop[5];
pr4[idx] = cell_prop[7];
}
}
}
}
}
LOGINFO("time for calculating PANPHASIA for level %d : %f s, %f µs/cell", level, get_wtime() - t1,
1e6 * (get_wtime() - t1) / ((double)nxremap[2] * (double)nxremap[1] * (double)nxremap[0]));
LOGUSER("\033[31mtiming level %d [adv_panphasia_cell_properties2]: %f s \033[0m", level, get_wtime() - tp);
tp = get_wtime();
/////////////////////////////////////////////////////////////////////////
// transform and convolve with Legendres
forward_transform_field(pr1, nxremap[0], nxremap[1], nxremap[2]);
forward_transform_field(pr2, nxremap[0], nxremap[1], nxremap[2]);
forward_transform_field(pr3, nxremap[0], nxremap[1], nxremap[2]);
forward_transform_field(pr4, nxremap[0], nxremap[1], nxremap[2]);
#pragma omp parallel for
for (int i = 0; i < nxremap[0]; i++)
for (int j = 0; j < nxremap[1]; j++)
for (int k = 0; k < nxremap[2] / 2 + 1; k++)
{
size_t idx = ((size_t)i * nxremap[1] + (size_t)j) * (nxremap[2] / 2 + 1) + (size_t)k;
double fx(1.0), fy(1.0), fz(1.0), arg = 0.;
complex gx(0., 0.), gy(0., 0.), gz(0., 0.);
int ii(i), jj(j), kk(k);
if (i > nxremap[0] / 2)
ii -= nxremap[0];
if (j > nxremap[1] / 2)
jj -= nxremap[1];
// int kkmax = std::max(abs(ii),std::max(abs(jj),abs(kk)));
if (ii != 0)
{
arg = M_PI * (double)ii / (double)nxremap[0];
fx = sin(arg) / arg;
gx = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
}
else
{
fx = 1.0;
gx = 0.0;
}
if (jj != 0)
{
arg = M_PI * (double)jj / (double)nxremap[1];
fy = sin(arg) / arg;
gy = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
}
else
{
fy = 1.0;
gy = 0.0;
}
if (kk != 0)
{
arg = M_PI * (double)kk / (double)nxremap[2];
fz = sin(arg) / arg;
gz = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
}
else
{
fz = 1.0;
gz = 0.0;
}
if (abs(ii) != nxremap[0] / 2 && abs(jj) != nxremap[1] / 2 &&
abs(kk) != nxremap[2] / 2)
{ // kkmax != nxremap[2]/2 ){
complex x, y1(RE(pc1[idx]), IM(pc1[idx])), y2(RE(pc2[idx]), IM(pc2[idx])), y3(RE(pc3[idx]), IM(pc3[idx])),
y4(RE(pc4[idx]), IM(pc4[idx]));
x = 3.0 * (y1 * gx * gy * fz + y2 * fx * gy * gz + y3 * gx * fy * gz) + sqrt(27.0) * y4 * gx * gy * gz;
RE(pc0[idx]) = RE(pc0[idx]) + x.real();
IM(pc0[idx]) = IM(pc0[idx]) + x.imag();
}
}
LOGUSER("\033[31mtiming level %d [build panphasia field2]: %f s\033[0m", level, get_wtime() - tp);
tp = get_wtime();
// END
//***************************************************************
// Compute Panphasia values of p011, p101, p110, p111 coefficients
// and combine with p000, p001, p010, p100 and indep field.
//****************************************************************
/////////////////////////////////////////////////////////////////////////
// do we need to cut off the small scales?
// int nn = 1<<level;
if (incongruent_fields_)
{
LOGUSER("Remapping fields from dimension %d -> %d", nxremap[0], nx_m[0]);
memset(pr1, 0, ngp * sizeof(fftw_real));
#pragma omp parallel for
for (int i = 0; i < nxremap[0]; i++)
for (int j = 0; j < nxremap[1]; j++)
for (int k = 0; k < nxremap[2] / 2 + 1; k++)
{
int ii = (i > nxremap[0] / 2) ? i - nxremap[0] : i, jj = (j > nxremap[1] / 2) ? j - nxremap[1] : j, kk = k;
int ia(abs(ii)), ja(abs(jj)), ka(abs(kk));
if (ia < nx_m[0] / 2 && ja < nx_m[1] / 2 && ka < nx_m[2] / 2)
{
size_t idx = ((size_t)(i)*nxremap[1] + (size_t)(j)) * (nxremap[2] / 2 + 1) + (size_t)(k);
int ir = (ii < 0) ? ii + nx_m[0] : ii, jr = (jj < 0) ? jj + nx_m[1] : jj, kr = kk; // never negative
size_t idx2 = ((size_t)ir * nx_m[1] + (size_t)jr) * ((size_t)nx_m[2] / 2 + 1) + (size_t)kr;
complex x(RE(pc0[idx]), IM(pc0[idx]));
double total_phase_shift;
total_phase_shift = inter_grid_phase_adjustment_ * (double)(ii + jj + kk);
x = x * exp(complex(0.0, total_phase_shift));
RE(pc1[idx2]) = x.real();
IM(pc1[idx2]) = x.imag();
}
}
memcpy(pr0, pr1, ngp * sizeof(fftw_real));
}
// if (level == 9)
// {
// LOGUSER("DC mode of level is %g", RE(pc0[0]));
// // RE(pc0[0]) = 1e8;
// // IM(pc0[0]) = 0.0;
// }
LOGUSER("\033[31mtiming level %d [remap noncongruent]: %f s\033[0m", level, get_wtime() - tp);
tp = get_wtime();
/////////////////////////////////////////////////////////////////////////
// transform back
backward_transform_field(pr0, nx_m[0], nx_m[1], nx_m[2]);
/////////////////////////////////////////////////////////////////////////
// copy to random data structure
delete[] pr1;
delete[] pr2;
delete[] pr3;
delete[] pr4;
LOGUSER("Copying random field data %d,%d,%d -> %d,%d,%d", nxremap[0], nxremap[1], nxremap[2], nx[0], nx[1], nx[2]);
// n = 1<<level;
// ng = n;
// ngp = ng*ng*2*(ng/2+1);
double sum = 0.0, sum2 = 0.0;
size_t count = 0;
/*double norm = 1.0 / sqrt((double)nxremap[0] * (double)nxremap[1] * (double)nxremap[2] * (double)nx[0] *
(double)nx[1] * (double)nx[2]);*/
double norm = 1.0 / sqrt((double)nxremap[0] * (double)nxremap[1] * (double)nxremap[2] * (double)nx_m[0] *
(double)nx_m[1] * (double)nx_m[2]);
#pragma omp parallel for reduction(+ \
: sum, sum2, count)
for (int k = 0; k < nx[2]; ++k) // ARJ - swapped roles of i,k, and reverse ordered loops
for (int j = 0; j < nx[1]; ++j)
for (int i = 0; i < nx[0]; ++i)
{
size_t idx = ((size_t)(i + iexpand_left[0]) * nx_m[1] + (size_t)(j + iexpand_left[1])) * (nx_m[2] + 2) + (size_t)(k + iexpand_left[2]);
R(i, j, k) = pr0[idx] * norm;
sum += R(i, j, k);
sum2 += R(i, j, k) * R(i, j, k);
++count;
}
delete[] pr0;
sum /= (double)count;
sum2 /= (double)count;
sum2 = (sum2 - sum * sum);
LOGUSER("done with PANPHASIA for level %d:\n mean=%g, var=%g", level, sum, sum2);
LOGUSER("Copying into R array: nx[0],nx[1],nx[2] %d %d %d \n", nx[0], nx[1], nx[2]);
LOGINFO("PANPHASIA level %d mean and variance are\n <p> = %g | var(p) = %g", level, sum, sum2);
}
namespace
{
RNG_plugin_creator_concrete<RNG_panphasia> creator("PANPHASIA");
}
#endif

10
src/poisson.cc
View file

@ -724,12 +724,16 @@ double fft_poisson_plugin::gradient(int dir, grid_hierarchy &u, grid_hierarchy &
}
#endif
/*double ktot = sqrt(ii*ii+jj*jj+k*k);
if( ktot >= nx/2 )//dir == 0 && i==nx/2 || dir == 1 && j==ny/2 || dir == 2 && k==nz/2 )
if( (dir == 0 && i==nx/2) || (dir == 1 && j==ny/2) || (dir == 2 && k==nz/2) )
{
#ifdef FFTW3
cdata[idx][0] = 0.0;
cdata[idx][1] = 0.0;
#else
cdata[idx].re = 0.0;
cdata[idx].im = 0.0;
}*/
#endif
}
}
RE(cdata[0]) = 0.0;

1831
src/random.cc
View file

File diff suppressed because it is too large Load diff

369
src/random.hh
View file

@ -4,7 +4,7 @@
a code to generate multi-scale initial conditions
for cosmological simulations
Copyright (C) 2010 Oliver Hahn
Copyright (C) 2010-23 by Oliver Hahn
*/
@ -13,15 +13,17 @@
// #define DEGRADE_RAND2
//.....................................
#ifndef __RANDOM_HH
#define __RANDOM_HH
#pragma once
#define DEF_RAN_CUBE_SIZE 32
#include <fstream>
#include <algorithm>
#include <map>
#ifdef _OPENMP
#include <omp.h>
#endif
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
@ -30,18 +32,23 @@
#include "mesh.hh"
#include "mg_operators.hh"
#include "constraints.hh"
// #include "convolution_kernel.hh"
#include "density_grid.hh"
class RNG_plugin
{
protected:
config_file *pcf_; //!< pointer to config_file from which to read parameters
const refinement_hierarchy *prefh_; //!< pointer to refinement hierarchy structure containing the grid sizes
public:
explicit RNG_plugin(config_file &cf)
: pcf_(&cf)
explicit RNG_plugin(config_file &cf) //, const refinement_hierarchy& refh )
: pcf_(&cf) //, prefh_( & refh )
{
}
virtual ~RNG_plugin() {}
virtual bool is_multiscale() const = 0;
virtual void fill_grid(int level, DensityGrid<real_t> &R) = 0;
virtual void initialize_for_grid_structure(const refinement_hierarchy &refh) = 0;
};
struct RNG_plugin_creator
@ -65,356 +72,54 @@ struct RNG_plugin_creator_concrete : public RNG_plugin_creator
}
//! create an instance of the plugin
RNG_plugin *create(config_file &cf) const
RNG_plugin *create(config_file &cf) const //, const refinement_hierarchy& refh ) const
{
return new Derived(cf);
return new Derived(cf); //, refh );
}
};
typedef RNG_plugin RNG_instance;
RNG_plugin *select_RNG_plugin(config_file &cf);
/*!
* @brief encapsulates all things random number generator related
*/
template <typename T>
class random_numbers
{
public:
unsigned
res_, //!< resolution of the full mesh
cubesize_, //!< size of one independent random number cube
ncubes_; //!< number of random number cubes to cover the full mesh
long baseseed_; //!< base seed from which cube seeds are computed
protected:
//! vector of 3D meshes (the random number cubes) with random numbers
std::vector<Meshvar<T> *> rnums_;
//! map of 3D indices to cube index
std::map<size_t, size_t> cubemap_;
typedef std::map<size_t, size_t>::iterator cubemap_iterator;
protected:
//! register a cube with the hash map
void register_cube(int i, int j, int k);
//! fills a subcube with random numbers
double fill_cube(int i, int j, int k);
//! subtract a constant from an entire cube
void subtract_from_cube(int i, int j, int k, double val);
//! copy random numbers from a cube to a full grid array
template <class C>
void copy_cube(int i, int j, int k, C &dat)
{
int offi, offj, offk;
offi = i * cubesize_;
offj = j * cubesize_;
offk = k * cubesize_;
i = (i + ncubes_) % ncubes_;
j = (j + ncubes_) % ncubes_;
k = (k + ncubes_) % ncubes_;
size_t icube = (i * ncubes_ + j) * ncubes_ + k;
cubemap_iterator it = cubemap_.find(icube);
if (it == cubemap_.end())
{
LOGERR("attempting to copy data from non-existing RND cube %d,%d,%d", i, j, k);
throw std::runtime_error("attempting to copy data from non-existing RND cube");
}
size_t cubeidx = it->second;
for (int ii = 0; ii < (int)cubesize_; ++ii)
for (int jj = 0; jj < (int)cubesize_; ++jj)
for (int kk = 0; kk < (int)cubesize_; ++kk)
dat(offi + ii, offj + jj, offk + kk) = (*rnums_[cubeidx])(ii, jj, kk);
}
//! free the memory associated with a subcube
void free_cube(int i, int j, int k);
//! initialize member variables and allocate memory
void initialize(void);
//! fill a cubic subvolume of the full grid with random numbers
double fill_subvolume(int *i0, int *n);
//! fill an entire grid with random numbers
double fill_all(void);
//! fill an external array instead of the internal field
template <class C>
double fill_all(C &dat)
{
double sum = 0.0;
for (int i = 0; i < (int)ncubes_; ++i)
for (int j = 0; j < (int)ncubes_; ++j)
for (int k = 0; k < (int)ncubes_; ++k)
{
int ii(i), jj(j), kk(k);
register_cube(ii, jj, kk);
}
#pragma omp parallel for reduction(+ \
: sum)
for (int i = 0; i < (int)ncubes_; ++i)
for (int j = 0; j < (int)ncubes_; ++j)
for (int k = 0; k < (int)ncubes_; ++k)
{
int ii(i), jj(j), kk(k);
ii = (ii + ncubes_) % ncubes_;
jj = (jj + ncubes_) % ncubes_;
kk = (kk + ncubes_) % ncubes_;
sum += fill_cube(ii, jj, kk);
copy_cube(ii, jj, kk, dat);
free_cube(ii, jj, kk);
}
return sum / (ncubes_ * ncubes_ * ncubes_);
}
//! write the number of allocated random number cubes to stdout
void print_allocated(void);
public:
//! constructor
random_numbers(unsigned res, unsigned cubesize, long baseseed, int *x0, int *lx);
//! constructor for constrained fine field
random_numbers(random_numbers<T> &rc, unsigned cubesize, long baseseed,
bool kspace = false, bool isolated = false, int *x0_ = NULL, int *lx_ = NULL, bool zeromean = true);
//! constructor
random_numbers(unsigned res, unsigned cubesize, long baseseed, bool zeromean = true);
//! constructor to read white noise from file
random_numbers(unsigned res, std::string randfname, bool rndsign);
//! copy constructor for averaged field (not copying) hence explicit!
explicit random_numbers(/*const*/ random_numbers<T> &rc, bool kdegrade = true);
//! destructor
~random_numbers()
{
for (unsigned i = 0; i < rnums_.size(); ++i)
if (rnums_[i] != NULL)
delete rnums_[i];
rnums_.clear();
}
//! access a random number, this allocates a cube and fills it with consistent random numbers
inline T &operator()(int i, int j, int k, bool fillrand = true)
{
int ic, jc, kc, is, js, ks;
if (ncubes_ == 0)
throw std::runtime_error("random_numbers: internal error, not properly initialized");
//... determine cube
ic = (int)((double)i / cubesize_ + ncubes_) % ncubes_;
jc = (int)((double)j / cubesize_ + ncubes_) % ncubes_;
kc = (int)((double)k / cubesize_ + ncubes_) % ncubes_;
size_t icube = ((size_t)ic * ncubes_ + (size_t)jc) * ncubes_ + (size_t)kc;
cubemap_iterator it = cubemap_.find(icube);
if (it == cubemap_.end())
{
LOGERR("Attempting to copy data from non-existing RND cube %d,%d,%d @ %d,%d,%d", ic, jc, kc, i, j, k);
throw std::runtime_error("attempting to copy data from non-existing RND cube");
}
size_t cubeidx = it->second;
if (rnums_[cubeidx] == NULL)
{
LOGERR("Attempting to access data from non-allocated RND cube %d,%d,%d", ic, jc, kc);
throw std::runtime_error("attempting to access data from non-allocated RND cube");
}
//... determine cell in cube
is = (i - ic * cubesize_ + cubesize_) % cubesize_;
js = (j - jc * cubesize_ + cubesize_) % cubesize_;
ks = (k - kc * cubesize_ + cubesize_) % cubesize_;
return (*rnums_[cubeidx])(is, js, ks);
}
//! free all cubes
void free_all_mem(void)
{
for (unsigned i = 0; i < rnums_.size(); ++i)
if (rnums_[i] != NULL)
{
delete rnums_[i];
rnums_[i] = NULL;
}
}
};
RNG_plugin *select_RNG_plugin(config_file &cf); //, const refinement_hierarchy& refh );
/*!
* @brief encapsulates all things for multi-scale white noise generation
*/
template <typename rng, typename T>
template <typename T>
class random_number_generator
{
protected:
config_file *pcf_;
refinement_hierarchy *prefh_;
constraint_set constraints;
int levelmin_,
levelmax_,
levelmin_seed_;
std::vector<long> rngseeds_;
std::vector<std::string> rngfnames_;
bool disk_cached_;
bool restart_;
std::vector<std::vector<T> *> mem_cache_;
unsigned ran_cube_size_;
protected:
//! checks if the specified string is numeric
bool is_number(const std::string &s);
//! parses the random number parameters in the conf file
void parse_rand_parameters(void);
//! correct coarse grid averages for the change in small scale when using Fourier interpolation
void correct_avg(int icoarse, int ifine);
//! the main driver routine for multi-scale white noise generation
void compute_random_numbers(void);
//! store the white noise fields in memory or on disk
void store_rnd(int ilevel, rng *prng);
// const refinement_hierarchy * prefh_;
RNG_plugin *generator_;
int levelmin_, levelmax_;
public:
//! constructor
random_number_generator(config_file &cf, refinement_hierarchy &refh, transfer_function *ptf = NULL);
random_number_generator(config_file &cf, transfer_function *ptf = NULL)
: pcf_(&cf) //, prefh_( &refh )
{
levelmin_ = pcf_->getValue<int>("setup", "levelmin");
levelmax_ = pcf_->getValue<int>("setup", "levelmax");
generator_ = select_RNG_plugin(cf);
}
//! destructor
~random_number_generator();
~random_number_generator()
{
}
//! initialize_for_grid_structure
void initialize_for_grid_structure(const refinement_hierarchy &refh)
{
generator_->initialize_for_grid_structure(refh);
}
//! load random numbers to a new array
template <typename array>
void load(array &A, int ilevel)
{
if (restart_)
LOGINFO("Attempting to restart using random numbers for level %d\n from file \'wnoise_%04d.bin\'.", ilevel, ilevel);
if (disk_cached_)
{
char fname[128];
sprintf(fname, "wnoise_%04d.bin", ilevel);
LOGUSER("Loading white noise from file \'%s\'...", fname);
std::ifstream ifs(fname, std::ios::binary);
if (!ifs.good())
{
LOGERR("White noise file \'%s\'was not found.", fname);
throw std::runtime_error("A white noise file was not found. This is an internal inconsistency and bad.");
}
int nx, ny, nz;
ifs.read(reinterpret_cast<char *>(&nx), sizeof(int));
ifs.read(reinterpret_cast<char *>(&ny), sizeof(int));
ifs.read(reinterpret_cast<char *>(&nz), sizeof(int));
if (nx != (int)A.size(0) || ny != (int)A.size(1) || nz != (int)A.size(2))
{
if (nx == (int)A.size(0) * 2 && ny == (int)A.size(1) * 2 && nz == (int)A.size(2) * 2)
{
std::cerr << "CHECKPOINT" << std::endl;
int ox = nx / 4, oy = ny / 4, oz = nz / 4;
std::vector<T> slice(ny * nz, 0.0);
for (int i = 0; i < nx; ++i)
{
ifs.read(reinterpret_cast<char *>(&slice[0]), ny * nz * sizeof(T));
if (i < ox)
continue;
if (i >= 3 * ox)
break;
#pragma omp parallel for
for (int j = oy; j < 3 * oy; ++j)
for (int k = oz; k < 3 * oz; ++k)
A(i - ox, j - oy, k - oz) = slice[j * nz + k];
}
ifs.close();
}
else
{
LOGERR("White noise file is not aligned with array. File: [%d,%d,%d]. Mem: [%d,%d,%d].",
nx, ny, nz, A.size(0), A.size(1), A.size(2));
throw std::runtime_error("White noise file is not aligned with array. This is an internal inconsistency and bad.");
}
}
else
{
for (int i = 0; i < nx; ++i)
{
std::vector<T> slice(ny * nz, 0.0);
ifs.read(reinterpret_cast<char *>(&slice[0]), ny * nz * sizeof(T));
#pragma omp parallel for
for (int j = 0; j < ny; ++j)
for (int k = 0; k < nz; ++k)
A(i, j, k) = slice[j * nz + k];
}
ifs.close();
}
}
else
{
LOGUSER("Copying white noise from memory cache...");
if (mem_cache_[ilevel - levelmin_] == NULL)
LOGERR("Tried to access mem-cached random numbers for level %d. But these are not available!\n", ilevel);
int nx(A.size(0)), ny(A.size(1)), nz(A.size(2));
if ((size_t)nx * (size_t)ny * (size_t)nz != mem_cache_[ilevel - levelmin_]->size())
{
LOGERR("White noise file is not aligned with array. File: [%d,%d,%d]. Mem: [%d,%d,%d].", nx, ny, nz, A.size(0), A.size(1), A.size(2));
throw std::runtime_error("White noise file is not aligned with array. This is an internal inconsistency and bad");
}
#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)
A(i, j, k) = (*mem_cache_[ilevel - levelmin_])[((size_t)i * ny + (size_t)j) * nz + (size_t)k];
std::vector<T>().swap(*mem_cache_[ilevel - levelmin_]);
delete mem_cache_[ilevel - levelmin_];
mem_cache_[ilevel - levelmin_] = NULL;
}
generator_->fill_grid(ilevel, A);
}
};
typedef random_numbers<real_t> rand_nums;
typedef random_number_generator<rand_nums, real_t> rand_gen;
#endif //__RANDOM_HH
using noise_generator = random_number_generator<real_t>;