mirror of
https://github.com/cosmo-sims/MUSIC.git
synced 2024-09-20 18:13:45 +02:00
WIP added panphasia files. does not work yet
This commit is contained in:
Oliver Hahn
parent
4aeb06b1d5
commit
040324f346
5 changed files with 4858 additions and 204 deletions
|
|
@ -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
683
ext/panphasia/generic_lecuyer.f90
Normal 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
3334
ext/panphasia/panphasia_routines.f
Normal file
File diff suppressed because it is too large
Load diff
|
|
@ -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
|
||||
821
src/plugins/random_panphasia.cc
Normal file
821
src/plugins/random_panphasia.cc
Normal file
|
|
@ -0,0 +1,821 @@
|
|||
#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_;
|
||||
|
||||
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());
|
||||
}
|
||||
|
||||
void initialize_for_grid_structure(const refinement_hierarchy &refh) {
|
||||
prefh_ = &refh;
|
||||
levelmin_ = prefh_->levelmin();
|
||||
levelmin_final_ = pcf_->getValue<unsigned>("setup", "levelmin");
|
||||
levelmax_ = prefh_->levelmax();
|
||||
|
||||
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{
|
||||
ileft_corner[k] = (ileft[k] - nx[k]/4 + (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];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
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]));
|
||||
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]));
|
||||
|
||||
//////////////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
LOGINFO("\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
|
||||
|
||||
LOGINFO("\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);
|
||||
|
||||
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();
|
||||
|
||||
// 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]));
|
||||
|
||||
LOGINFO("\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();
|
||||
}
|
||||
}
|
||||
|
||||
LOGINFO("\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_) {
|
||||
|
||||
LOGINFO("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 ){
|
||||
LOGINFO("DC mode of level is %g",RE(pc0[0]));
|
||||
//RE(pc0[0]) = 1e8;
|
||||
//IM(pc0[0]) = 0.0;
|
||||
}
|
||||
|
||||
LOGINFO("\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;
|
||||
|
||||
LOGINFO("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
|
||||
Loading…
Reference in a new issue