diff --git a/CMakeLists.txt b/CMakeLists.txt index e83719d..51ed9bb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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}) diff --git a/ext/panphasia/generic_lecuyer.f90 b/ext/panphasia/generic_lecuyer.f90 new file mode 100644 index 0000000..13f53ed --- /dev/null +++ b/ext/panphasia/generic_lecuyer.f90 @@ -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 +! + + +!}}} + diff --git a/ext/panphasia/panphasia_routines.f b/ext/panphasia/panphasia_routines.f new file mode 100644 index 0000000..2e1bfbd --- /dev/null +++ b/ext/panphasia/panphasia_routines.f @@ -0,0 +1,3334 @@ +c=====================================================================================c +c +c The code below was written by: Adrian Jenkins, +c Institute for Computational Cosmology +c Department of Physics +c South Road +c Durham, DH1 3LE +c United Kingdom +c +c This file is part of the software made public in +c Jenkins and Booth 2013 - arXiv:1306.XXXX +c +c The software computes the Panphasia Gaussian white noise field +c realisation described in detail in Jenkins 2013 - arXiv:1306.XXXX +c +c +c +c This software is free, subject to a agreeing licence conditions: +c +c +c (i) you will publish the phase descriptors and reference Jenkins (13) +c for any new simulations that use Panphasia phases. You will pass on this +c condition to others for any software or data you make available publically +c or privately that makes use of Panphasia. +c +c (ii) that you will ensure any publications using results derived from Panphasia +c will be submitted as a final version to arXiv prior to or coincident with +c publication in a journal. +c +c (iii) that you report any bugs in this software as soon as confirmed to +c A.R.Jenkins@durham.ac.uk +c +c (iv) that you understand that this software comes with no warranty and that is +c your responsibility to ensure that it is suitable for the purpose that +c you intend. +c +c=====================================================================================c + +c===================================================================================== +c List of subroutines and arguments. Each of these is documented in c +c arXiV/1306.XXXX c +c c +c Adrian Jenkins, 24/6/2013. c +c------------------------------------------------------------------------------------- +c Version 1.000 +c=================================================================================== + + module pan_state + use Rand + implicit none + integer maxdim_, maxlev_, maxpow_ + parameter (maxdim_=60,maxlev_=50, maxpow_ = 3*maxdim_) + integer nmulti_ + parameter (nmulti_=64) + integer range_max + parameter(range_max=10000) + integer indmin,indmax + parameter (indmin=-1, indmax=60) + + + type state_data + integer base_state(5), base_lev_start(5,0:maxdim_) + TYPE(Rand_offset) :: poweroffset(0:maxpow_) + TYPE(Rand_offset) :: superjump + TYPE(Rand_state) :: current_state(-1:maxpow_) + + integer layer_min,layer_max,indep_field + +! This module stores information needed to access the part of Panphasia +! selected by a particular descriptor. + integer*8 xorigin_store(0:1,0:1,0:1) + integer*8 yorigin_store(0:1,0:1,0:1) + integer*8 zorigin_store(0:1,0:1,0:1) + + integer*4 lev_common + integer*4 layer_min_store,layer_max_store + + integer*8 ix_abs_store,iy_abs_store,iz_abs_store + integer*8 ix_per_store,iy_per_store,iz_per_store + integer*8 ix_rel_store,iy_rel_store,iz_rel_store + + real*8 exp_coeffs(8,0:7,-1:maxdim_) + integer*8 xcursor(0:maxdim_),ycursor(0:maxdim_),zcursor(0:maxdim_) + +c Local box parameters + + integer*4 ixshift(0:1,0:1,0:1) + integer*4 iyshift(0:1,0:1,0:1) + integer*4 izshift(0:1,0:1,0:1) + + +c more state variables + real*8 cell_data(9,0:7) + integer*4 ixh_last,iyh_last,izh_last + integer init + + integer return_cell_props_init + integer reset_lecuyer_state_init + integer*8 p_xcursor(indmin:indmax),p_ycursor(indmin:indmax),p_zcursor(indmin:indmax) + + + + end type state_data + + + +c Switch for enabling custom spherical function +c Set isub_spherical_function = 1 to turn on the spherical function + integer*4 isub_spherical_function + parameter (isub_spherical_function=0) + + end module pan_state + + +c================================================================================ +c Begin white noise routines +c================================================================================ + recursive subroutine start_panphasia(ldata,descriptor,ngrid,VERBOSE) + use pan_state + implicit none + type(state_data), intent(inout) :: ldata + character*100 descriptor + integer ngrid + integer VERBOSE + + + + integer*4 wn_level_base,i_base,i_base_y,i_base_z + integer*8 i_xorigin_base,i_yorigin_base,i_zorigin_base, check_rand + character*20 name + + integer ratio + integer lextra + integer level_p + + + integer*8 ix_abs,iy_abs,iz_abs + integer*8 ix_per,iy_per,iz_per + integer*8 ix_rel,iy_rel,iz_rel + + !integer layer_min,layer_max,indep_field + !common /oct_range/ layer_min,layer_max,indep_field + + call parse_descriptor(descriptor ,wn_level_base,i_xorigin_base,i_yorigin_base, + & i_zorigin_base,i_base,i_base_y,i_base_z,check_rand,name) + + + lextra = (log10(real(ngrid)/real(i_base))+0.001)/log10(2.0) + ratio = 2**lextra + + if (ratio*i_base.ne.ngrid) + &stop 'Value of ngrid inconsistent with dim of region in Panphasia' + + level_p = wn_level_base + lextra + + ix_abs = ishft(i_xorigin_base,lextra) + iy_abs = ishft(i_yorigin_base,lextra) + iz_abs = ishft(i_zorigin_base,lextra) + + ix_per = i_base*ratio + iy_per = i_base*ratio + iz_per = i_base*ratio + +c Set the refinement position at the origin. + + ix_rel = 0 + iy_rel = 0 + iz_rel = 0 + + call set_phases_and_rel_origin(ldata,descriptor,level_p,ix_rel,iy_rel,iz_rel,VERBOSE) + +c Finally set the octree functions required for making cosmological +c initial conditions. These are passed using a common block. + + ldata%layer_min = 0 + ldata%layer_max = level_p + ldata%indep_field = 1 + + end +c================================================================================= + recursive subroutine set_phases_and_rel_origin(ldata,descriptor,lev,ix_rel,iy_rel,iz_rel,VERBOSE) + use pan_state + !use descriptor_phases + implicit none + type(state_data), intent(inout) :: ldata + character*100 descriptor + integer lev + integer*8 ix_abs,iy_abs,iz_abs + integer*8 ix_per,iy_per,iz_per + integer*8 ix_rel,iy_rel,iz_rel + integer*8 xorigin,yorigin,zorigin + + integer VERBOSE + integer MYID + integer*8 maxco + integer i + integer px,py,pz + + integer lnblnk + integer*8 mconst + parameter(mconst = 2147483647_Dint) + + integer*4 wn_level_base,i_base,i_base_y,i_base_z + integer*8 i_xorigin_base,i_yorigin_base,i_zorigin_base, check_rand + integer lextra,ratio + character*20 phase_name + +c----------------------------------------------------------------------------------------------- + + call initialise_panphasia(ldata) + + call validate_descriptor(ldata, descriptor,-1,check_rand) + + call parse_descriptor(descriptor ,wn_level_base,i_xorigin_base,i_yorigin_base, + & i_zorigin_base,i_base,i_base_y,i_base_z,check_rand,phase_name) + lextra = lev - wn_level_base + ratio = 2**lextra + + ix_abs = ishft(i_xorigin_base,lextra) + iy_abs = ishft(i_yorigin_base,lextra) + iz_abs = ishft(i_zorigin_base,lextra) + + ix_per = i_base*ratio + iy_per = i_base*ratio + iz_per = i_base*ratio + +c------------------------------------------------------------------------- +c Error checking +c------------------------------------------------------------------------- + if ((lev.lt.0).or.(lev.gt.maxlev_)) stop 'Level out of range! (1)' + + + maxco = 2_dint**lev + + if (ix_abs.lt.0) stop 'Error: ix_abs negative (1)' + if (iy_abs.lt.0) stop 'Error: iy_abs negative (1)' + if (iz_abs.lt.0) stop 'Error: iz_abs negative (1)' + + if (ix_rel.lt.0) stop 'Error: ix_rel negative (1)' + if (iy_rel.lt.0) stop 'Error: iy_rel negative (1)' + if (iz_rel.lt.0) stop 'Error: iz_rel negative (1)' + + + if (ix_abs+ix_rel.ge.maxco) + & stop 'Error: ix_abs + ix_rel out of range. (1)' + if (iy_abs+iy_rel.ge.maxco) + & stop 'Error: iy_abs + iy_rel out of range. (1)' + if (iz_abs+iz_rel.ge.maxco) + & stop 'Error: iz_abs + iz_rel out of range. (1)' + +c---------------------------------------------------------------------------------------- +c To allow the local box to wrap around, if needed, define a series of eight +c 'origins'. For many purposes (ix,iy,iz) = (0,0,0) is the only origin needed. + + + do px=0,1 + do py=0,1 + do pz=0,1 + + xorigin = max(0,( ix_abs + ix_rel - px*ix_per )/2) + yorigin = max(0,( iy_abs + iy_rel - py*iy_per )/2) + zorigin = max(0,( iz_abs + iz_rel - pz*iz_per )/2) + + ldata%ixshift(px,py,pz) = max(0, ix_abs + ix_rel -px*ix_per) - 2*xorigin + ldata%iyshift(px,py,pz) = max(0, iy_abs + iy_rel -py*iy_per) - 2*yorigin + ldata%izshift(px,py,pz) = max(0, iz_abs + iz_rel -pz*iz_per) - 2*zorigin + + +c Store box details: store the positions at level lev-1 + + + ldata%xorigin_store(px,py,pz) = xorigin + ldata%yorigin_store(px,py,pz) = yorigin + ldata%zorigin_store(px,py,pz) = zorigin + + enddo + enddo + enddo + + ldata%lev_common = lev + + + ldata%ix_abs_store = ix_abs + ldata%iy_abs_store = iy_abs + ldata%iz_abs_store = iz_abs + + ldata%ix_per_store = ix_per + ldata%iy_per_store = iy_per + ldata%iz_per_store = iz_per + + ldata%ix_rel_store = ix_rel + ldata%iy_rel_store = iy_rel + ldata%iz_rel_store = iz_rel + + +c Reset all cursor values to negative numbers. + + do i=0,maxdim_ + ldata%xcursor(i) = -999 + ldata%ycursor(i) = -999 + ldata%zcursor(i) = -999 + enddo + if (VERBOSE.gt.1) then + if (MYID.lt.1) then + print*,'----------------------------------------------------------' + print*,'Successfully initialised Panphasia box at level ',lev + write (6,105) ix_abs,iy_abs,iz_abs + write (6,106) ix_rel,iy_rel,iz_rel + write (6,107) ix_per,iy_per,iz_per + write (6,*) 'Phases used: ',descriptor(1:lnblnk(descriptor)) + print*,'----------------------------------------------------------' + endif + endif + 105 format(' Abs origin: (',i12,',',i12,',',i12,')') + 106 format(' Rel origin: (',i12,',',i12,',',i12,')') + 107 format(' Periods : (',i12,',',i12,',',i12,')') + end +c================================================================================ + recursive subroutine initialise_panphasia( ldata ) + use Rand + use pan_state + implicit none + + type(state_data), intent(inout) :: ldata + + TYPE(Rand_state) :: state + TYPE(Rand_offset) :: offset + integer ninitialise + parameter (ninitialise=218) + integer i + real*8 rand_num + + + call Rand_seed(state,ninitialise) + + call Rand_save(ldata%base_state,state) + + call Rand_set_offset(offset,1) + +c Calculate offsets of powers of 2 times nmulti +c + + do i=0,maxpow_ + ldata%poweroffset(i) = Rand_mul_offset(offset,nmulti_) + offset = Rand_mul_offset(offset,2) + enddo + + +c Compute the base state for each level. + + call Rand_load(state,ldata%base_state) + state = Rand_step(state,8) + + do i=0,maxdim_ + call Rand_save(ldata%base_lev_start(1,i),state) + state = Rand_boost(state,ldata%poweroffset(3*i)) + enddo + +c Set superjump to value 2**137 - used occasionally in computing Gaussian variables +c when the value of the returned random number is less an 10-6. + + call Rand_set_offset(ldata%superjump,1) + + do i=1,137 + ldata%superjump = Rand_mul_offset(ldata%superjump,2) + enddo + + +c Run time test to see if one particular value can be recovered. + + call Rand_load(state,ldata%base_lev_start(1,34)) + call Rand_real(rand_num,state) + + if (abs(rand_num- 0.828481889948473d0).gt.1.e-14) then + print*,'Error in initialisation!' + print*,'Rand_num = ',rand_num + print*,'Target value = ', 0.828481889948473d0 + stop + endif + return + end +c================================================================================= + recursive subroutine panphasia_cell_properties(ldata,ixcell,iycell,izcell,cell_prop) + use pan_state + implicit none + type(state_data), intent(inout) :: ldata + !integer layer_min,layer_max,indep_field + !common /oct_range/ layer_min,layer_max,indep_field + integer*4 ixcell,iycell,izcell + real*8 cell_prop(9) + + call adv_panphasia_cell_properties(ldata,ixcell,iycell,izcell,ldata%layer_min, + & ldata%layer_max,ldata%indep_field,cell_prop) + return + end +c================================================================================= + recursive subroutine adv_panphasia_cell_properties(ldata,ixcell,iycell,izcell,layer_min, + & layer_max,indep_field,cell_prop) + use pan_state + !use descriptor_phases + implicit none + + type(state_data), intent(inout) :: ldata + + integer*4 lev + integer*4 ixcell,iycell,izcell + integer layer_min,layer_max,indep_field + real*8 cell_prop(9) +c real*8 cell_data(9,0:7) + integer*4 j,l,lx,ly,lz + integer*4 px,py,pz + +c integer*4 ixh_last,iyh_last,izh_last + +c integer init +c data init/0/ +c save init,cell_data,ixh_last,iyh_last,izh_last ! Keep internal state + + integer*4 ixh,iyh,izh + + lev = ldata%lev_common + +c------- Error checking ----------------------------- + + if (layer_min.gt.layer_max) then + + if (layer_min-layer_max.eq.1) then ! Not necessarily bad. No octree basis functions + do j=1,9 ! required at this level and position. + cell_prop(j) = 0.0d0 ! Set returned cell_prop data to zero. + enddo + return + endif + + print*,'Warning: layer_min.gt.layer_max!' + print*,'layer_min = ',layer_min + print*,'layer_max = ',layer_max + print*,'ixcell,iycell,izcell',ixcell,iycell,izcell + + call flush(6) + stop 'Error: layer_min.gt.layer_max' + endif + + if (layer_max.gt.ldata%lev_common) then + print*,'lev_common = ',ldata%lev_common + print*,'layer_min = ',layer_min + print*,'layer_max = ',layer_max + stop 'Error: layer_max.gt.lev_common' + endif + if ((indep_field.lt.-1).or.(indep_field.gt.1)) + & stop 'Error: indep_field out of range' + +c---------------------------------------------------- +c Check which 'origin' to use. + + px = 0 + py = 0 + pz = 0 + + if (ldata%ix_rel_store+ixcell.ge.ldata%ix_per_store) px = 1 ! Crossed x-periodic bndy + if (ldata%iy_rel_store+iycell.ge.ldata%iy_per_store) py = 1 ! Crossed y-periodic bndy + if (ldata%iz_rel_store+izcell.ge.ldata%iz_per_store) pz = 1 ! Crossed z-periodic bndy +c---------------------------------------------------- + + + ixh = (ixcell+ldata%ixshift(px,py,pz) )/2 + iyh = (iycell+ldata%iyshift(px,py,pz) )/2 + izh = (izcell+ldata%izshift(px,py,pz) )/2 + + lx = mod(ixcell+ldata%ixshift(px,py,pz) ,2) + ly = mod(iycell+ldata%iyshift(px,py,pz) ,2) + lz = mod(izcell+ldata%izshift(px,py,pz) ,2) + + + l = 4*lx + 2*ly + lz ! Determine which cell is required + +cc------------------ If no new evalation is needed skip assignment ----- + if ((ldata%init.eq.1).and.(ixh.eq.ldata%ixh_last).and.(iyh.eq.ldata%iyh_last).and. + & (izh.eq.ldata%izh_last).and.(layer_min.eq.ldata%layer_min_store).and. + & (layer_max.eq.ldata%layer_max_store)) goto 24 +cc----------------------------------------------------------------------------- + + + call return_cell_props(ldata,lev,ixh,iyh,izh,px,py,pz,layer_min, + & layer_max,indep_field,ldata%cell_data) + +c Remember previous values. + + ldata%ixh_last = ixh + ldata%iyh_last = iyh + ldata%izh_last = izh + + + 24 continue + + + do j=1,9 + cell_prop(j) = ldata%cell_data(j,l) ! Copy the required data + enddo + + if (ldata%init.eq.0) ldata%init=1 + + return + end +c================================================================================= + recursive subroutine return_cell_props(ldata,lev_input,ix_half,iy_half,iz_half, + & px,py,pz,layer_min,layer_max,indep_field,cell_data) + use Rand + use pan_state + !use descriptor_phases + implicit none + type(state_data), intent(inout) :: ldata + integer lev_input,ix_half,iy_half,iz_half,px,py,pz + integer layer_min,layer_max,indep_field + real*8 cell_data(9,0:7) + + real*8 garray(0:63) + integer lev + integer*8 xarray,yarray,zarray + + integer i,istart,icell_name + + +c integer init +c data init/0/ +c save init + + + +c-------------------------------------------------------- +c--------------------------- Initialise level -1 -------- +c-------------------------------------------------------- + + if (ldata%return_cell_props_init.eq.0) then ! First time called. Set up the Legendre coefficients + ldata%return_cell_props_init = 1 ! for the root cell. This is the first term on the + call Rand_load(ldata%current_state(-1),ldata%base_state) ! right hand side of the equation in appendix C of + call return_gaussian_array(ldata,-1,8,garray) ! Jenkins 2013 that defines PANPHASIA. + ldata%exp_coeffs(1,0,-1) = garray(0) + ldata%exp_coeffs(2,0,-1) = garray(1) + ldata%exp_coeffs(3,0,-1) = garray(2) + ldata%exp_coeffs(4,0,-1) = garray(3) + ldata%exp_coeffs(5,0,-1) = garray(4) + ldata%exp_coeffs(6,0,-1) = garray(5) + ldata%exp_coeffs(7,0,-1) = garray(6) + ldata%exp_coeffs(8,0,-1) = garray(7) + + ldata%layer_min_store = layer_min + ldata%layer_max_store = layer_max + + endif + +c-------------------------------------------------------- +c---------------------------- Error checking ------------ +c-------------------------------------------------------- + + lev = lev_input-1 + + if (lev_input.ne.ldata%lev_common) stop 'Box initialised at a different level !' + if (ix_half.lt.0) then + print*,'ix_half negative',ix_half + stop 'ix_half out of range!' + endif + if (iy_half.lt.0) stop 'iy_half out of range!' + if (iz_half.lt.0) then + print*,'iz_half negative',iz_half + stop 'iz_half out of range!' + endif + + + xarray = ldata%xorigin_store(px,py,pz) + ix_half + yarray = ldata%yorigin_store(px,py,pz) + iy_half + zarray = ldata%zorigin_store(px,py,pz) + iz_half + + +c If layer_max or layer_min have changed, rebuild from the start and reset the +c recorded value of layer_max and layer_min + + if ((layer_max.ne.ldata%layer_max_store).or.(layer_min.ne.ldata%layer_min_store)) then + + if (layer_min.gt.layer_max) stop 'layer_min > layer_max : 2' + + istart = max(1,layer_min-1) + + ldata%layer_max_store = layer_max + ldata%layer_min_store = layer_min + + goto 10 + + endif + + + if ((xarray.eq.ldata%xcursor(lev)).and.(yarray.eq.ldata%ycursor(lev)).and.(zarray.eq.ldata%zcursor(lev))) return ! Nothing to do. + +c=========================================================================================================== +c------------- First determine which levels need to be (re)computed +c=========================================================================================================== + + istart = 0 + do i=lev-1,0,-1 + if ((ishft(xarray,i-lev).eq.ldata%xcursor(i)).and.(ishft(yarray,i-lev).eq.ldata%ycursor(i)).and. + & (ishft(zarray,i-lev).eq.ldata%zcursor(i))) then + istart = i+1 + goto 10 + endif + enddo + + 10 continue + + +c==================================================================================== +c------------- Now compute each level as required and update (x,y,z) cursor variables +c==================================================================================== + + do i=istart,lev + + icell_name = 0 + + ldata%xcursor(i) = ishft(xarray,i-lev) + ldata%ycursor(i) = ishft(yarray,i-lev) + ldata%zcursor(i) = ishft(zarray,i-lev) + + if (btest(ldata%xcursor(i),0)) icell_name = icell_name + 4 + if (btest(ldata%ycursor(i),0)) icell_name = icell_name + 2 + if (btest(ldata%zcursor(i),0)) icell_name = icell_name + 1 + + call reset_lecuyer_state(ldata,i,ldata%xcursor(i),ldata%ycursor(i),ldata%zcursor(i)) + + if (isub_spherical_function.ne.1) then + call return_gaussian_array(ldata,i,64,garray) + else + call return_oct_sf_expansion(ldata,i,lev,ldata%xcursor(i),ldata%ycursor(i),ldata%zcursor(i), + & 64,garray) + endif + + + call evaluate_panphasia(ldata,i,maxdim_,garray,layer_min, + & layer_max, indep_field, icell_name,cell_data,ldata%exp_coeffs) + + enddo + return + end +c================================================================================= + recursive subroutine evaluate_panphasia(ldata,nlev,maxdim,g, + & layer_min,layer_max,indep_field,icell_name,cell_data,leg_coeff) + use pan_state + implicit none +c--------------------------------------------------------------------------------- +c This subroutine calculates the Legendre block coefficients for the eight child +c cells of an octree cell. +c +c----------------- Define subroutine arguments ----------------------------------- + type(state_data), intent(inout) :: ldata + integer nlev,maxdim + integer layer_min,layer_max,indep_field + integer icell_name + real*8 leg_coeff(0:7,0:7,-1:maxdim),cell_data(0:8,0:7) + real*8 g(*) + +c----------------- Define constants using notation from appendix A of Jenkins 2013 + + real*8 a1,a2,b1,b2,b3,c1,c2,c3,c4 + + parameter(a1 = 0.5d0*sqrt(3.0d0), a2 = 0.5d0) + + parameter(b1 = 0.75d0, b2 = 0.25d0*sqrt(3.0d0)) + parameter(b3 = 0.25d0) + + parameter(c1 = sqrt(27.0d0/64.0d0), c2 = 0.375d0) + parameter(c3 = sqrt(3.0d0/64.0d0), c4 = 0.125d0) + +c----------------- Define octree variables -------------------------------- + + real*8 coeff_p000, coeff_p001, coeff_p010, coeff_p011 + real*8 coeff_p100, coeff_p101, coeff_p110, coeff_p111 + + real*8 positive_octant_lc(0:7,0:1,0:1,0:1),temp_value(0:7,0:7) + integer i,j,ix,iy,iz + integer icx,icy,icz + integer iox,ioy,ioz + real*8 parity,isig + real*8 usually_rooteighth_factor +c-------------------------------------------------------------------------- + +c------------- Set the Legendre block coefficients for the parent cell +c itself. These are either inherited from the octree above +c or set to zero depending on which levels of the octree +c have been selected to be populated with the octree +c basis functions. +c--------------------------------------------------------------------------- + if (nlev.ge.layer_min) then + coeff_p000 = leg_coeff(0,icell_name,nlev-1) + coeff_p001 = leg_coeff(1,icell_name,nlev-1) + coeff_p010 = leg_coeff(2,icell_name,nlev-1) + coeff_p011 = leg_coeff(3,icell_name,nlev-1) + coeff_p100 = leg_coeff(4,icell_name,nlev-1) + coeff_p101 = leg_coeff(5,icell_name,nlev-1) + coeff_p110 = leg_coeff(6,icell_name,nlev-1) + coeff_p111 = leg_coeff(7,icell_name,nlev-1) + else + coeff_p000 = 0.0d0 + coeff_p001 = 0.0d0 + coeff_p010 = 0.0d0 + coeff_p011 = 0.0d0 + coeff_p100 = 0.0d0 + coeff_p101 = 0.0d0 + coeff_p110 = 0.0d0 + coeff_p111 = 0.0d0 + endif + +c Apply layer_max and indep_field inputs --------------------------------- + + if (indep_field.ne.-1) then + usually_rooteighth_factor = sqrt(0.125d0) + else + usually_rooteighth_factor = 0.0d0 ! This option returns only the indep field. + endif ! For use in testing only. + + if (nlev.ge.layer_max) then + do i=1,56 + g(i) = 0.0d0 ! Set octree coefficients to zero as not required. + enddo + endif + + if (indep_field.eq.0) then ! Set the independent field to zero as not required. + do i=57,64 + g(i) = 0.0d0 + enddo + endif +c----------------------------------------------------------------------------- +c +c +c The calculations immediately below evalute the eight Legendre block coefficients for the +c child cell that is furthest from the absolute coordiate origin of the octree - we call +c this the positive octant cell. +c +c The coefficients are given by a set of matrix equations which combine the +c coefficients of the Legendre basis functions of the parent cell itself, with +c the coefficients from the octree basis functions that occupy the +c parent cell. +c +c The Legendre basis function coefficients of the parent cell are stored in +c the variables, coeff_p000 - coeff_p111 and are initialise above. +c +c The coefficients of the octree basis functions are determined by the +c first 56 entries of the array g, which is passed down into this +c subroutine. +c +c These two sources of information are combined using a set of linear equations. +c The coefficients of these linear equations are taken from the inverses or +c equivalently transposes of the matrices given in appendix A of Jenkins 2013. +c The matrices in appendix A define the PANPHASIA octree basis functions +c in terms of Legendre blocks. +c +c All of the Legendre block functions of the parent cell, and the octree basis +c functions of the parent cell share one of eight distinct symmetries with respect to +c reflection about the x1=0,x2=0,x3=0 planes (where the origin is taken as the parent +c cell centre and x1,x2,x3 are parallel to the cell edges). +c +c Each function has either purely reflectional symmetry (even parity) or +c reflectional symmetry with a sign change (odd parity) about each of the three principal +c planes through the cell centre. There are therefore 8 parity types. We can label each +c parity type with a binary triplet. So 000 is pure reflectional symmetry about +c all of the principal planes. +c +c In the code below the parent cell Legendre block functions, and octree functions are +c organised into eight groups each with eight members. Each group has a common +c parity type. +c +c We keep the contributions of each parity type to each of the eight Legendre basis +c functions occupying the positive octant cell separate. Once they have all been +c computed, we can apply the different symmetry operations and determine the +c Legendre block basis functions for all eight child cells at the same time. +c--------------------------------------------------------------------------------------- +c 000 parity + + positive_octant_lc(0, 0,0,0) = 1.0d0*coeff_p000 + positive_octant_lc(1, 0,0,0) = -1.0d0*g(1) + positive_octant_lc(2, 0,0,0) = -1.0d0*g(2) + positive_octant_lc(3, 0,0,0) = 1.0d0*g(3) + positive_octant_lc(4, 0,0,0) = -1.0d0*g(4) + positive_octant_lc(5, 0,0,0) = 1.0d0*g(5) + positive_octant_lc(6, 0,0,0) = 1.0d0*g(6) + positive_octant_lc(7, 0,0,0) = -1.0d0*g(7) + +c 100 parity + + positive_octant_lc(0, 1,0,0) = a1*coeff_p100 - a2*g(8) + positive_octant_lc(1, 1,0,0) = g(9) + positive_octant_lc(2, 1,0,0) = g(10) + positive_octant_lc(3, 1,0,0) = -g(11) + positive_octant_lc(4, 1,0,0) = a2*coeff_p100 + a1*g(8) + positive_octant_lc(5, 1,0,0) = -g(12) + positive_octant_lc(6, 1,0,0) = -g(13) + positive_octant_lc(7, 1,0,0) = g(14) + +c 010 parity + + positive_octant_lc(0, 0,1,0) = a1*coeff_p010 - a2*g(15) + positive_octant_lc(1, 0,1,0) = g(16) + positive_octant_lc(2, 0,1,0) = a2*coeff_p010 + a1*g(15) + positive_octant_lc(3, 0,1,0) = -g(17) + positive_octant_lc(4, 0,1,0) = g(18) + positive_octant_lc(5, 0,1,0) = -g(19) + positive_octant_lc(6, 0,1,0) = -g(20) + positive_octant_lc(7, 0,1,0) = g(21) + + +c 001 parity + + positive_octant_lc(0, 0,0,1) = a1*coeff_p001 - a2*g(22) + positive_octant_lc(1, 0,0,1) = a2*coeff_p001 + a1*g(22) + positive_octant_lc(2, 0,0,1) = g(23) + positive_octant_lc(3, 0,0,1) = -g(24) + positive_octant_lc(4, 0,0,1) = g(25) + positive_octant_lc(5, 0,0,1) = -g(26) + positive_octant_lc(6, 0,0,1) = -g(27) + positive_octant_lc(7, 0,0,1) = g(28) + +c 110 parity + + positive_octant_lc(0, 1,1,0) = b1*coeff_p110 - b2*g(29) + b3*g(30) - b2*g(31) + positive_octant_lc(1, 1,1,0) = -g(32) + positive_octant_lc(2, 1,1,0) = b2*coeff_p110 - b3*g(29) - b2*g(30) + b1*g(31) + positive_octant_lc(3, 1,1,0) = g(33) + positive_octant_lc(4, 1,1,0) = b2*coeff_p110 + b1*g(29) + b2*g(30) + b3*g(31) + positive_octant_lc(5, 1,1,0) = g(34) + positive_octant_lc(6, 1,1,0) = b3*coeff_p110 + b2*g(29) - b1*g(30) - b2*g(31) + positive_octant_lc(7, 1,1,0) = -g(35) + + +c 011 parity + + positive_octant_lc(0, 0,1,1) = b1*coeff_p011 - b2*g(36) + b3*g(37) - b2*g(38) + positive_octant_lc(1, 0,1,1) = b2*coeff_p011 - b3*g(36) - b2*g(37) + b1*g(38) + positive_octant_lc(2, 0,1,1) = b2*coeff_p011 + b1*g(36) + b2*g(37) + b3*g(38) + positive_octant_lc(3, 0,1,1) = b3*coeff_p011 + b2*g(36) - b1*g(37) - b2*g(38) + positive_octant_lc(4, 0,1,1) = -g(39) + positive_octant_lc(5, 0,1,1) = g(40) + positive_octant_lc(6, 0,1,1) = g(41) + positive_octant_lc(7, 0,1,1) = -g(42) + +c 101 parity + + positive_octant_lc(0, 1,0,1) = b1*coeff_p101 - b2*g(43) + b3*g(44) - b2*g(45) + positive_octant_lc(1, 1,0,1) = b2*coeff_p101 - b3*g(43) - b2*g(44) + b1*g(45) + positive_octant_lc(2, 1,0,1) = -g(46) + positive_octant_lc(3, 1,0,1) = g(47) + positive_octant_lc(4, 1,0,1) = b2*coeff_p101 + b1*g(43) + b2*g(44) + b3*g(45) + positive_octant_lc(5, 1,0,1) = b3*coeff_p101 + b2*g(43) - b1*g(44) - b2*g(45) + positive_octant_lc(6, 1,0,1) = g(48) + positive_octant_lc(7, 1,0,1) = -g(49) + +c 111 parity + + positive_octant_lc(0, 1,1,1) = c1*coeff_p111 - c2*g(50) - c2*g(51) - c2*g(52) + c3*g(53) + c3*g(54) + c3*g(55) - c4*g(56) + positive_octant_lc(1, 1,1,1) = c2*coeff_p111 + c1*g(50) - c2*g(51) + c2*g(52) - c3*g(53) + c3*g(54) + c4*g(55) + c3*g(56) + positive_octant_lc(2, 1,1,1) = c2*coeff_p111 + c2*g(50) + c1*g(51) - c2*g(52) - c3*g(53) - c4*g(54) + c3*g(55) - c3*g(56) + positive_octant_lc(3, 1,1,1) = c3*coeff_p111 - c3*g(50) - c3*g(51) + c4*g(52) - c1*g(53) - c2*g(54) - c2*g(55) - c2*g(56) + positive_octant_lc(4, 1,1,1) = c2*coeff_p111 - c2*g(50) + c2*g(51) + c1*g(52) + c4*g(53) - c3*g(54) + c3*g(55) + c3*g(56) + positive_octant_lc(5, 1,1,1) = c3*coeff_p111 + c3*g(50) - c4*g(51) - c3*g(52) + c2*g(53) - c1*g(54) - c2*g(55) + c2*g(56) + positive_octant_lc(6, 1,1,1) = c3*coeff_p111 + c4*g(50) + c3*g(51) + c3*g(52) + c2*g(53) + c2*g(54) - c1*g(55) - c2*g(56) + positive_octant_lc(7, 1,1,1) = c4*coeff_p111 - c3*g(50) + c3*g(51) - c3*g(52) - c2*g(53) + c2*g(54) - c2*g(55) + c1*g(56) +c-------------------------------------------------------------------------------------------- +c +c +c We now calculate the Legendre basis coefficients for all eight child cells +c by applying the appropriate reflectional parities to the coefficients +c calculated above for the positive octant child cell. +c +c See equations A2 and A3 in appendix A of Jenkins 2013. +c +c The reflectional parity is given by (ix,iy,iz) loops below. +c +c The (icx,icy,icz) loops below, loop over the eight child cells. +c +c The positive octant child cell is given below by (icx=icy=icz=0) or i=7. +c +c The combination ix*icx +iy*icy +iz*icz is either even or odd, depending +c on whether the parity change is even or odd. +c +c The variables iox,ioy,ioz are used to loop over the different +c types of Legendre basis function. +c +c The combination iox*icx + ioy*icy + ioz*icz is either even and odd +c and identifies which coefficients keep or change sign respectively +c due to a pure reflection about the principal planes. +c-------------------------------------------------------------------------------------------- + + do iz=0,7 + do iy=0,7 + temp_value(iy,iz) = 0.0d0 ! Zero temporary sums + enddo + enddo +c-------------------------------------------------------------------------------------------- + do iz=0,1 ! Loop over z parity (0=keep sign, 1=change sign) + do iy=0,1 ! Loop over y parity (0=keep sign, 1=change sign) + do ix=0,1 ! Loop over x parity (0=keep sign, 1=change sign) + + + do icx=0,1 ! Loop over x-child cells + do icy=0,1 ! Loop over y-child cells + do icz=0,1 ! Loop over z-child cells + + if (mod(ix*icx+iy*icy+iz*icz,2).eq.0) then + parity = 1.0d0 + else + parity =-1.0d0 + endif + + i = 7 - 4*icx -2*icy - icz ! Calculate which child cell this is. + + + do iox=0,1 ! Loop over Legendre basis function type + do ioy=0,1 ! Loop over Legendre basis function type + do ioz=0,1 ! Loop over Legendre basis function type + + j = 4*iox + 2*ioy + ioz + + if (mod(iox*icx + ioy*icy + ioz*icz,2).eq.0) then + isig = parity + else + isig = -parity + endif + + temp_value(j,i) = temp_value(j,i) + isig*positive_octant_lc(j,ix,iy,iz) + + enddo + enddo + enddo + + enddo + enddo + enddo + + enddo + enddo + enddo + + +c Assign values of the output variables + + do i=0,7 + do j=0,7 + leg_coeff(j,i,nlev) = temp_value(j,i)*usually_rooteighth_factor + cell_data(j,i) = leg_coeff(j,i,nlev) + enddo + enddo + +c Finally set the independent field values + + cell_data(8,0) = g(57) + cell_data(8,1) = g(58) + cell_data(8,2) = g(59) + cell_data(8,3) = g(60) + cell_data(8,4) = g(61) + cell_data(8,5) = g(62) + cell_data(8,6) = g(63) + cell_data(8,7) = g(64) + + + return + end +c================================================================================= + recursive subroutine reset_lecuyer_state(ldata,lev,xcursor,ycursor,zcursor) + use pan_state + implicit none + + type(state_data), intent(inout) :: ldata + integer lev + integer*8 xcursor,ycursor,zcursor + +c integer indmin,indmax +c parameter (indmin=-1, indmax=60) +c integer*8 p_xcursor(indmin:indmax),p_ycursor(indmin:indmax),p_zcursor(indmin:indmax) +c save p_xcursor,p_ycursor,p_zcursor + integer i +c integer init +c data init/0/ +c save init + + if (ldata%reset_lecuyer_state_init.eq.0) then ! Initialise p_cursor variables with + ldata%reset_lecuyer_state_init = 1 ! negative values. + do i=indmin,indmax + ldata%p_xcursor(i) = -9999 + ldata%p_ycursor(i) = -9999 + ldata%p_zcursor(i) = -9999 + enddo + endif + + if ( (xcursor.eq.ldata%p_xcursor(lev)).and.(ycursor.eq.ldata%p_ycursor(lev)).and. + & (zcursor.eq.ldata%p_zcursor(lev)+1)) then + ldata%p_xcursor(lev) = xcursor + ldata%p_ycursor(lev) = ycursor + ldata%p_zcursor(lev) = zcursor + return + endif + + call advance_current_state(ldata,lev,xcursor,ycursor,zcursor) + + ldata%p_xcursor(lev) = xcursor + ldata%p_ycursor(lev) = ycursor + ldata%p_zcursor(lev) = zcursor + + + return + end +c================================================================================= + recursive subroutine advance_current_state(ldata,lev,x,y,z) + use Rand + use pan_state + !use descriptor_phases + implicit none + + type(state_data), intent(inout) :: ldata + + integer lev + integer*8 x,y,z + + integer*8 lev_range + + TYPE(Rand_offset) :: offset1,offset2 + TYPE(Rand_offset) :: offset_x,offset_y,offset_z,offset_total + + integer ndiv,nrem + integer*8 ndiv8,nrem8 + integer nfactor + parameter (nfactor=291071) ! Value unimportant except has to be > 262144 + + +c----- First some error checking ------------------------------------------ + if ((lev.lt.0).or.(lev.gt.maxlev_)) stop 'Level out of range! (2)' + + lev_range = 2_dint**lev + + + if ((x.lt.0).or.(x.ge.lev_range)) then + print*,'x,lev,lev_range',x,lev,lev_range + call flush(6) + stop 'x out of range!' + endif + if ((y.lt.0).or.(y.ge.lev_range)) then + print*,'y,lev,lev_range',y,lev,lev_range + stop 'y out of range!' + endif + if ((z.lt.0).or.(z.ge.lev_range)) stop 'z out of range!' +c---------------------------------------------------------------------------- +c +c Note the Rand_set_offset subroutine takes an integer*4 value +c for the offset value. For this reason we need to use integer*4 +c values - ndiv,nrem. As a precaution an explicit check is made +c to be sure that these values are calculated correctly. +c--------------------------------------------------------------------------- + + + call Rand_load(ldata%current_state(lev),ldata%base_lev_start(1,lev)) + + if (lev.eq.0) return + +c Calculate z-offset + + ndiv = z/nfactor + nrem = z - ndiv*nfactor + ndiv8 = ndiv + nrem8 = nrem + + if (ndiv8*nfactor+nrem8.ne.z) stop 'Error in z ndiv nrem' + + call Rand_set_offset(offset1,ndiv) + offset1 = Rand_mul_offset(offset1,nfactor) + call Rand_set_offset(offset2,nrem) + offset2 = Rand_add_offset(offset1,offset2) + offset_z = Rand_mul_offset(offset2,nmulti_) + +c Calculate y-offset + + ndiv = y/nfactor + nrem = y - ndiv*nfactor + ndiv8 = ndiv + nrem8 = nrem + + if (ndiv8*nfactor+nrem8.ne.y) stop 'Error in y ndiv nrem' + + offset1 = Rand_mul_offset(ldata%poweroffset(lev),ndiv) + offset1 = Rand_mul_offset(offset1,nfactor) + offset2 = Rand_mul_offset(ldata%poweroffset(lev),nrem) + offset_y = Rand_add_offset(offset1,offset2) + +c Calculate x-offset + + ndiv = x/nfactor + nrem = x - ndiv*nfactor + ndiv8 = ndiv + nrem8 = nrem + + if (ndiv8*nfactor+nrem8.ne.x) then + print*,'ndiv,nfactor,nrem,x',ndiv,nfactor,nrem,x + print*,'ndiv*nfactor+nrem',ndiv*nfactor+nrem + print*,'x-ndiv*nfactor-nrem',x-ndiv*nfactor-nrem + stop 'Error in x ndiv nrem' + endif + + offset1 = Rand_mul_offset(ldata%poweroffset(2*lev),ndiv) + offset1 = Rand_mul_offset(offset1,nfactor) + offset2 = Rand_mul_offset(ldata%poweroffset(2*lev),nrem) + offset_x = Rand_add_offset(offset1,offset2) + + offset1 = Rand_add_offset(offset_x,offset_y) + offset_total = Rand_add_offset(offset1, offset_z) + + ldata%current_state(lev) = Rand_boost(ldata%current_state(lev),offset_total) + + return + end +c================================================================================= + recursive subroutine return_gaussian_array(ldata,lev,ngauss,garray) + use Rand + use pan_state + implicit none + type(state_data), intent(inout) :: ldata + integer lev,ngauss + real*8 garray(0:*) + TYPE(Rand_state) :: state + real*8 PI + parameter (PI=3.1415926535897932384d0) + real*8 branch + parameter (branch=1.d-6) + integer iloop + + real*8 temp,mag,ang + integer i + + if (mod(ngauss,2).ne.0) + & stop 'Error in return_gaussian_array - even pairs only' + +c First obtain a set of uniformly distributed pseudorandom numbers +c between 0 and 1. The method used is described in detail in +c appendix B of Jenkins 2013. + + do i=0,ngauss-1 + call Rand_real(garray(i),ldata%current_state(lev)) + + if (garray(i).lt.branch) then + garray(i) = branch + state = Rand_boost(ldata%current_state(lev),ldata%superjump) + iloop = 0 + 10 continue + call Rand_real(temp,state) + iloop = iloop+1 + if (temp.lt.branch) then + garray(i) = garray(i)*branch + state = Rand_boost(state,ldata%superjump) + if (iloop.gt.100) then + print*,'Too may iterations in return_gaussian_array!' + call flush(6) + stop + endif + goto 10 + else + garray(i) = garray(i)*temp + endif + endif + enddo + +c Apply Box-Muller transformation to create pairs of Gaussian +c pseudorandom numbers. + + do i=0,ngauss/2-1 + + mag = sqrt(-2.0d0*log(garray(2*i))) + ang = 2.0d0*PI*garray(2*i+1) + + garray(2*i) = mag*cos(ang) + garray(2*i+1) = mag*sin(ang) + + enddo + end +c================================================================================= + recursive subroutine parse_descriptor(string,l,ix,iy,iz,side1,side2,side3,check_int,name) + implicit none + integer nchar + parameter(nchar=100) + character*100 string + integer*4 l,side1,side2,side3,ierror + integer*8 ix,iy,iz + integer*8 check_int + character*20 name + + + integer i,ip,iq,ir + + ierror = 0 + + ip = 1 + do while (string(ip:ip).eq.' ') + ip = ip + 1 + enddo + + if (string(ip:ip+7).ne.'[Panph1,') then + ierror = 1 + print*,string(ip:ip+7) + goto 10 + endif + + ip = ip+8 + if (string(ip:ip).ne.'L') then + ierror = 2 + goto 10 + endif + + ip = ip+1 + + iq = ip + scan( string(ip:nchar),',') -1 + + if (ip.eq.iq) then + ierror = 3 + goto 10 + endif + + + read (string(ip:iq),*) l + + ip = iq+1 + + if (string(ip:ip).ne.'(') then + ierror = 4 + goto 10 + endif + + ip = ip+1 + + iq = ip + scan( string(ip:nchar),')') -2 + + read(string(ip:iq),*) ix,iy,iz + + ip = iq+2 + + if (string(ip:ip).ne.',') then + ierror = 5 + goto 10 + endif + + ip = ip+1 + if ((string(ip:ip).ne.'S').and.(string(ip:ip).ne.'D')) then + ierror = 6 + goto 10 + endif + + if (string(ip:ip).eq.'S') then + ip = ip + 1 + iq = ip + scan( string(ip:nchar),',') -2 + read (string(ip:iq),*) side1 + side2 = side1 + side3 = side1 + iq = iq+1 + if (string(iq:iq+2).ne.',CH') then + print*,string(ip:iq),string(iq:iq+2) + ierror = 6 + goto 10 + endif + else + ip = ip + 1 + if (string(ip:ip).ne.'(') then + ierror = 7 + goto 10 + endif + + + ip = ip + 1 + iq = ip + scan( string(ip:nchar),')') -2 + read (string(ip:iq),*) side1,side2,side3 + + iq = iq + 1 + + if (string(iq:iq).ne.')') then + ierror = 8 + goto 10 + endif + + iq = iq + 1 + + if (string(iq:iq+2).ne.',CH') then + ierror = 9 + goto 10 + endif + + endif + + ip = iq + 3 + + iq = ip + scan( string(ip:nchar),',') -2 + + read (string(ip:iq),*) check_int + + ip = iq + 1 + + if (string(ip:ip).ne.',') then + ierror = 10 + goto 10 + endif + + ip = ip+1 + + ir = ip + scan( string(ip:nchar),']') -2 + + iq = min(ir,ip+19) + + do i=1,20 + name(i:i)=' ' + enddo + + do i=ip,iq + name(i-ip+1:i-ip+1) = string(i:i) + enddo + + iq = ir + 1 + + if (string(iq:iq).ne.']') then + ierror = 11 + goto 10 + endif + + + 10 continue + + if (ierror.eq.0) return + + print*,'Error reading panphasian descriptor. Error number:',ierror + stop + + return + end +c================================================================================= + recursive subroutine compose_descriptor(l,ix,iy,iz,side,check_int,name,string) + implicit none + integer nchar + parameter(nchar=100) + character*100,intent(out)::string + character*20 name + integer*4 l,ltemp + integer*8 side + integer*8 ix,iy,iz + integer*8 check_int + + character*50 temp1,temp2,temp3,temp4,temp5,temp6 + integer lnblnk + + integer ip1,ip2,ip3,ip4,ip5,ip6 + + ltemp = l + + 5 continue + if ((mod(ix,2).eq.0).and.(mod(iy,2).eq.0).and.(mod(iz,2).eq.0).and.(mod(side,2).eq.0)) then + ix = ix/2 + iy = iy/2 + iz = iz/2 + side = side/2 + ltemp = ltemp-1 + goto 5 + endif + + + write (temp1,*) ltemp + ip1= scan(temp1,'0123456789') + write (temp2,*) ix + ip2= scan(temp2,'0123456789') + write (temp3,*) iy + ip3= scan(temp3,'0123456789') + write (temp4,*) iz + ip4= scan(temp4,'0123456789') + write (temp5,*) side + ip5= scan(temp5,'0123456789') + write (temp6,*) check_int + ip6= scan(temp6,'-0123456789') + + + string='[Panph1,L'//temp1(ip1:lnblnk(temp1))//',('//temp2(ip2:lnblnk(temp2)) + & //','//temp3(ip3:lnblnk(temp3))//','//temp4(ip4:lnblnk(temp4))//'),S' + & // temp5(ip5:lnblnk(temp5))//',CH'//temp6(ip6:lnblnk(temp6))// + & ','//name(1:lnblnk(name))//']' + + return + + end +c================================================================================= + recursive subroutine validate_descriptor(ldata,string,MYID,check_number) + use pan_state + implicit none + + type(state_data), intent(inout) :: ldata + character*100 string + integer*8 check_number + integer MYID + + character*20 phase_name + integer*4 lev + + integer*8 ix_abs,iy_abs,iz_abs + integer*4 ix_base,iy_base,iz_base + + + integer*8 xval,yval,zval + integer val_state(5) + + TYPE(Rand_state) :: state + + real*8 rand_num + integer*8 mconst,check_total,check_rand + parameter(mconst = 2147483647_Dint) + integer ascii_list(0:255) + integer*8 maxco + integer i + integer*8 ii + integer lnblnk + + + + call parse_descriptor(string,lev,ix_abs,iy_abs,iz_abs, + & ix_base,iy_base,iz_base,check_rand,phase_name) + +c------------------------------------------------------------------------- +c Some basic checking +c------------------------------------------------------------------------- + if ((lev.lt.0).or.(lev.gt.maxlev_)) then + print*,'lev,maxlev',lev,maxlev_ + call flush(6) + stop 'Level out of range! (3)' + endif + + if ((mod(ix_abs,2).eq.0).and.(mod(iy_abs,2).eq.0).and.(mod(iz_abs,2).eq.0).and. + & (mod(ix_base,2).eq.0).and.(mod(iy_base,2).eq.0).and.(mod(iz_base,2).eq.0)) + & stop 'Parameters not at lowest level' + + + maxco = 2_dint**lev + + if (ix_abs.lt.0) stop 'Error: ix_abs negative (2)' + if (iy_abs.lt.0) stop 'Error: iy_abs negative (2)' + if (iz_abs.lt.0) stop 'Error: iz_abs negative (2)' + + + if (ix_abs+ix_base.ge.maxco) + & stop 'Error: ix_abs + ix_per out of range.' + if (iy_abs+iy_base.ge.maxco) + & stop 'Error: iy_abs + iy_per out of range.' + if (iz_abs+iz_base.ge.maxco) + & stop 'Error: iz_abs + iz_per out of range.' + + check_total = 0 + + call initialise_panphasia(ldata) +c First corner + xval = ix_abs + ix_base - 1 + yval = iy_abs + zval = iz_abs + call advance_current_state(ldata,lev,xval,yval,zval) + call Rand_real(rand_num,ldata%current_state(lev)) + call Rand_save(val_state,ldata%current_state(lev)) + check_total = check_total + val_state(5) + if (MYID.eq.0) print*,'--------------------------------------' + if (MYID.eq.0) print*,'X-corner rand = ',rand_num + if (MYID.eq.0) print*,'State:',val_state +c Second corner + xval = ix_abs + yval = iy_abs + iy_base - 1 + zval = iz_abs + call advance_current_state(ldata,lev,xval,yval,zval) + call Rand_real(rand_num,ldata%current_state(lev)) + call Rand_save(val_state,ldata%current_state(lev)) + check_total = check_total + val_state(5) + if (MYID.eq.0) print*,'Y-corner rand = ',rand_num + if (MYID.eq.0) print*,'State:',val_state +c Third corner + xval = ix_abs + yval = iy_abs + zval = iz_abs + iz_base - 1 + call advance_current_state(ldata,lev,xval,yval,zval) + call Rand_real(rand_num,ldata%current_state(lev)) + call Rand_save(val_state,ldata%current_state(lev)) + check_total = check_total + val_state(5) + if (MYID.eq.0) print*,'z-corner rand = ',rand_num + if (MYID.eq.0) print*,'State:',val_state + if (MYID.eq.0) print*,'--------------------------------------' + +c Now encode the name. An integer for each ascii character is generated +c starting from the state which gives r0 - the first random number in +c Panphasia. The integer is in the range 0 - m-1. +c After making the list, then loop over non-blank characters +c in the name and take the ascii value, and sum the associated numbers. +c To avoid simple anagrams giving the same score, weight the integer +c by position in the string. Finally take mod m - to give the +c check number. + + call Rand_load(state,ldata%base_state) + + do i=0,255 + call Rand_real(rand_num,state) + call Rand_save(val_state,state) + ascii_list(i) = val_state(5) + enddo + + + + do ii=1,lnblnk(phase_name) + check_total = check_total + ii*ascii_list(iachar(phase_name(ii:ii))) + enddo + + + check_total = mod(check_total,mconst) + if (check_rand.eq.-999) then ! override the safety check number. + check_number = check_total + return + else + if (check_rand.ne.check_total) then + print*,'Inconsistency in the input panphasia descriptor ',MYID + print*,'Check_rand = ',check_rand + print*,'val_state(5) =',val_state(5) + print*,'xval,yval,zval',xval,yval,zval + print*,'lev_val = ',lev + call flush(6) + stop + endif + endif + + + return + end +c================================================================================= + recursive subroutine generate_random_descriptor(ldata,string) + use Rand + use pan_state + implicit none + type(state_data), intent(inout) :: ldata + character*100 string + character*100 instring + character*20 name + integer*4 unix_timestamp + + real*8 lbox + real*8 lpanphasia + parameter (lpanphasia = 25000000.0) ! Units of Mpc/h + integer level + integer*8 cell_dim + integer val_state(5) + + TYPE(Rand_state) :: state + TYPE(Rand_offset) :: offset + + real*8 rand_num1,rand_num2 + integer*8 mconst,check_int + parameter(mconst = 2147483647_Dint) + integer*8 mfac,imajor,iminor + parameter(mfac=33554332_Dint) + integer ascii_list(0:255) + integer i,lnblnk + integer*8 ii + integer mult + + integer*8 ixco,iyco,izco,irange + + print*,'___________________________________________________________' + print* + print*,' Generate a random descriptor ' + print* + print*,'The code uses the time (the unix timestamp) plus some extra ' + print*,'information personal to the user to choose a random region ' + print*,'within PANPHASIA. The user must also specify the side length' + print*,'of the cosmological volume. The code assumes that the whole of' + print*,'PANPHASIA is 25000 Gpc/h on a side and selects an appropriate ' + print*,'level in the octree for the descriptor. ' + print*,'Assuming this scaling the small scale power is defined down ' + print*,'to a mass scale of around 10^{-12} solar masses.' + print* + print*,'The user must also specify a human readable label for the ' + print*,'descriptor of less than 21 characters.' + print*,'___________________________________________________________' + print* + print*,'Press return to continue ' + read (*,*) + print* + print*,'___________________________________________________________' + print*,'Enter the box side-length in Mpc/h units' + read (*,*) lbox + print*,'___________________________________________________________' + print* + print* + 5 continue + print*,'Enter up to 20 character name to label the descriptor (no spaces)' + read (*,'(a)') name + if ((len_trim(instring).lt.21).or.(scan(name,' ').le.len_trim(name))) goto 5 + print*,'___________________________________________________________' + print* + print* + print*,'___________________________________________________________' + print*,'The phases for the simulation are described by whole octree ' + print*,'cells. Enter an odd integer that defines the number of cells ' + print*,'you require in one dimension. Choose this number carefully ' + print*,'as it will limit the possible 1-D sizes of the of the Fourier ' + print*,'transforms that can be used to make initial conditions to a product ' + print*,'of this integer times any power of two. In which case the only' + print*,'choice is 1.)' + print*,'(I would recommend 3 unless the initial condition code is' + print*,'incapable of using grid sizes that are not purely powers of two.' + print*,'___________________________________________________________' + print* + 7 continue + print*,'Enter number of octree cells on an edge (positive odd number only) ' + read (*,*) cell_dim + if ((cell_dim.le.0).or.(mod(cell_dim,2).eq.0)) goto 7 + print*,'___________________________________________________________' + call system('date +%s>tempfile_42526037646') + open(16,file='tempfile_42526037646',status='old') + read (16,*) unix_timestamp + close(16) + call system('/bin/rm tempfile_42526037646') + + print*,'Unix_timestamp determined. Value: ',unix_timestamp + print*,'___________________________________________________________' + print* + print* + print* + print*,'___________________________________________________________' + print*,'The code has just read the unix timestamp and will use this' + print*,'to help choose a random region in PANPHASIA. Although it is' + print*,'perhaps unlikely that someone else is also running this code at ' + print*,'the same time to the nearest second, to make it more likely' + print*,' still that the desciptor to be generated is unique' + print*,'please enter your name or some other piece of information' + print*,'below that you think is unlikely to be used by anyone else' + print*,'___________________________________________________________' + + print* + + 10 continue + print*,'Please enter your name (a minimum of six characters)' + read (*,'(a)') instring !' + if (len_trim(instring).lt.6) goto 10 + + level = int(log10(dble(cell_dim)*lpanphasia/lbox)/log10(2.0d0)) + + if (level.gt.50) stop 'level >50 ' + + + +c 'd' lines allow the generation of a large set of +c descriptors. Use to check that they are randomly +c positioned over the available volume. + + +c First use the unix timestamp to initialises the +c random generator. + + call Rand_seed(state,unix_timestamp) + + call Rand_save(ldata%base_state,state) + + +c First generate an integer from the user data. + call Rand_load(state,ldata%base_state) + + do i=0,255 + call Rand_real(rand_num1,state) + call Rand_save(val_state,state) + ascii_list(i) = val_state(5) + enddo + + call Rand_set_offset(offset,1) + + do ii=1,lnblnk(instring) + mult = mod(ii*ascii_list(iachar(instring(ii:ii))),mconst) + offset = Rand_mul_offset(offset,mult) + enddo + + call Rand_load(state,ldata%base_state) + state = Rand_boost(state,offset) ! Starting point for choosing location. + + 20 continue + + irange = 2_Dint**level + imajor = irange/mfac + iminor = mod(irange,mfac) + + call Rand_real(rand_num1,state) + call Rand_real(rand_num2,state) + + ixco = int(rand_num1*imajor)*mfac + int(rand_num2*iminor) + + if (ixco+cell_dim.ge.irange) goto 20 ! Invalid descriptor + + call Rand_real(rand_num1,state) + call Rand_real(rand_num2,state) + + iyco = int(rand_num1*imajor)*mfac + int(rand_num2*iminor) + + if (iyco+cell_dim.ge.irange) goto 20 ! Invalid descriptor + + call Rand_real(rand_num1,state) + call Rand_real(rand_num2,state) + + izco = int(rand_num1*imajor)*mfac + int(rand_num2*iminor) + + if (izco+cell_dim.ge.irange) goto 20 ! Invalid descriptor + + +c Value of the check digit is not known. Use validate_descriptor to compute it. + + check_int = -999 ! Special value required to make validate_descriptor + ! return the check digit. + + call compose_descriptor(level,ixco,iyco,izco,cell_dim,check_int,name,string) + + call validate_descriptor(ldata,string,-1,check_int) + + call compose_descriptor(level,ixco,iyco,izco,cell_dim,check_int,name,string) + + + return + end +c================================================================================= + recursive subroutine demo_basis_function_allocator + + implicit none + integer nmax + parameter (nmax=10) + + integer*4 wn_level(nmax) + + integer*8 ix_abs(nmax),iy_abs(nmax),iz_abs(nmax) + integer*8 ix_per(nmax),iy_per(nmax),iz_per(nmax) + integer*8 ix_rel(nmax),iy_rel(nmax),iz_rel(nmax) + integer*8 ix_dim(nmax),iy_dim(nmax),iz_dim(nmax) + + integer ix,iy,iz,nref + integer layer_min,layer_max,indep_field + + + integer*8 itot_int,itot_ib + + integer inv_open + +c Assign some trial values + + nref = 3 + inv_open=9 + + wn_level(1) = 22 + + ix_abs(1) = 2000000 + iy_abs(1) = 1500032 + iz_abs(1) = 2500032 + + ix_per(1) = 768 + iy_per(1) = 768 + iz_per(1) = 768 + + ix_rel(1) = 0 + iy_rel(1) = 0 + iz_rel(1) = 0 + + ix_dim(1) = 768 + iy_dim(1) = 768 + iz_dim(1) = 768 + + + wn_level(2) = 23 + + ix_abs(2) = 4000000 + iy_abs(2) = 3000064 + iz_abs(2) = 5000064 + + ix_per(2) = 1536 + iy_per(2) = 1536 + iz_per(2) = 1536 + + ix_rel(2) = 256 + iy_rel(2) = 16 + iz_rel(2) = 720 + + ix_dim(2) = 768 + iy_dim(2) = 768 + iz_dim(2) = 768 + + + wn_level(3) = 24 + + ix_abs(3) = 8000000 + iy_abs(3) = 6000128 + iz_abs(3) = 10000128 + + ix_per(3) = 3072 + iy_per(3) = 3072 + iz_per(3) = 3072 + + ix_rel(3) = 896 + iy_rel(3) = 432 + iz_rel(3) = 1840 + + ix_dim(3) = 768 + iy_dim(3) = 768 + iz_dim(3) = 768 + + + itot_int = 0 + itot_ib = 0 + + + + + open(10,file='ascii_dump_r1',status='unknown') + + ix=320 + do iy=0,767 + do iz=0,767 + call layer_choice(ix,iy,iz,1,nref,ix_abs,iy_abs,iz_abs, + & ix_per,iy_per,iz_per,ix_rel,iy_rel,iz_rel,ix_dim,iy_dim,iz_dim, + & wn_level,inv_open,layer_min,layer_max,indep_field) + write(10,*) iy,iz,layer_min,layer_max,indep_field + enddo + enddo + close(10) + + open(10,file='ascii_dump_r2',status='unknown') + + ix=384 + do iy=0,767 + do iz=0,767 + call layer_choice(ix,iy,iz,2,nref,ix_abs,iy_abs,iz_abs, + & ix_per,iy_per,iz_per,ix_rel,iy_rel,iz_rel,ix_dim,iy_dim,iz_dim, + & wn_level,inv_open,layer_min,layer_max,indep_field) + write(10,*) iy,iz,layer_min,layer_max,indep_field + enddo + enddo + close(10) + + open(10,file='ascii_dump_r3',status='unknown') + + ix=384 + do iy=0,767 + do iz=0,767 + call layer_choice(ix,iy,iz,3,nref,ix_abs,iy_abs,iz_abs, + & ix_per,iy_per,iz_per,ix_rel,iy_rel,iz_rel,ix_dim,iy_dim,iz_dim, + & wn_level,inv_open,layer_min,layer_max,indep_field) + write(10,*) iy,iz,layer_min,layer_max,indep_field + enddo + enddo + close(10) + end +c================================================================================= + recursive subroutine layer_choice(ix0,iy0,iz0,iref,nref, + & ix_abs,iy_abs,iz_abs,ix_per,iy_per,iz_per, + & ix_rel,iy_rel,iz_rel,ix_dim,iy_dim,iz_dim, + & wn_level,x_fact,layer_min,layer_max,indep_field) + implicit none + + integer ix0,iy0,iz0,iref,nref,isize,ibase + integer ix,iy,iz,irefplus + integer ione + + integer*8 ix_abs(nref),iy_abs(nref),iz_abs(nref) + integer*8 ix_per(nref),iy_per(nref),iz_per(nref) + integer*8 ix_rel(nref),iy_rel(nref),iz_rel(nref) + integer*8 ix_dim(nref),iy_dim(nref),iz_dim(nref) + + integer wn_level(nref) + integer layer_min,layer_max,indep_field,x_fact + integer idebug + + + integer interior,iboundary + + if (iref.eq.9999) then + idebug = 1 + else + idebug = 0 + endif + + ione = 1 + + irefplus = min(iref+1,nref) + + if (nref.eq.1) then ! Deal with simplest case + layer_min = 0 + layer_max = wn_level(1) + indep_field = 1 + if (idebug.eq.1) print*,'return 1' + return + endif + +c----------- Case of the top periodic refinement. For this refinement layer_min=0 as +c----------- all the larger basis functions must be included. By default layer_max +c----------- is set to wn_level(1) so all basis functions are included. A check is +c----------- made to determine if the lowest basis function can be included in the +c----------- next refinement. If it can the same process is repeated for the next +c----------- largest basis function and this is repeated until a failure occurs. + + if ((iref.eq.1).and.(nref.gt.1)) then + ibase = 1 + 10 continue + + ix = ishft(ishft(ix_abs(iref)+ix_rel(iref)+ix0,-ibase),ibase)-ix_abs(iref)-ix_rel(iref) + iy = ishft(ishft(iy_abs(iref)+iy_rel(iref)+iy0,-ibase),ibase)-iy_abs(iref)-iy_rel(iref) + iz = ishft(ishft(iz_abs(iref)+iz_rel(iref)+iz0,-ibase),ibase)-iz_abs(iref)-iz_rel(iref) + isize = ishft(ione,ibase) + + call inref(ix,iy,iz,isize,iref,irefplus,nref,wn_level, + & ix_abs,iy_abs,iz_abs,ix_per,iy_per,iz_per, + & ix_rel,iy_rel,iz_rel,ix_dim,iy_dim,iz_dim,x_fact, + & interior,iboundary) + + if ((interior.eq.1).and.(iboundary.eq.1)) then + ibase = ibase + 1 + goto 10 + endif + + layer_min = 0 + layer_max = wn_level(iref) - ibase + 1 + if (layer_max.ne.wn_level(iref)) then + indep_field = 0 + else + indep_field = 1 + endif + + if (idebug.eq.1) then + print*,'iref,wn_level(iref)',iref,wn_level(iref) + print*,'Return 2',layer_min,layer_max,indep_field + endif + + return + endif +c------------------------------------------------------------------------------------------ +c------------------------------------------------------------------------------------------ + + +c----------- For second or higher refinement determine layer_min by reference +c----------- to itself. In this case the loop continues until a basis function +c------------ is found which fits in a larger refinement + + ibase = 1 + + 20 continue + + + ix = ishft(ishft(ix_abs(iref)+ix_rel(iref)+ix0,-ibase),ibase)-ix_abs(iref)-ix_rel(iref) + iy = ishft(ishft(iy_abs(iref)+iy_rel(iref)+iy0,-ibase),ibase)-iy_abs(iref)-iy_rel(iref) + iz = ishft(ishft(iz_abs(iref)+iz_rel(iref)+iz0,-ibase),ibase)-iz_abs(iref)-iz_rel(iref) + isize = ishft(ione,ibase) + + call inref(ix,iy,iz,isize,iref,iref,nref,wn_level, + & ix_abs,iy_abs,iz_abs,ix_per,iy_per,iz_per, + & ix_rel,iy_rel,iz_rel,ix_dim,iy_dim,iz_dim,x_fact, + & interior,iboundary) + + if ((interior.eq.1).and.(iboundary.eq.1)) then + ibase = ibase + 1 + goto 20 + endif + + layer_min = wn_level(iref) - max(ibase-2,0) ! Take last suitable refinement + + +c----------- For an intermediate refinement define layer_max by reference to +c----------- the next refinement + + if (iref.lt.nref) then + ibase = 1 + + 30 continue + + ix = ishft(ishft(ix_abs(iref)+ix_rel(iref)+ix0,-ibase),ibase)-ix_abs(iref)-ix_rel(iref) + iy = ishft(ishft(iy_abs(iref)+iy_rel(iref)+iy0,-ibase),ibase)-iy_abs(iref)-iy_rel(iref) + iz = ishft(ishft(iz_abs(iref)+iz_rel(iref)+iz0,-ibase),ibase)-iz_abs(iref)-iz_rel(iref) + isize = ishft(ione,ibase) + + call inref(ix,iy,iz,isize,iref,irefplus,nref,wn_level, + & ix_abs,iy_abs,iz_abs,ix_per,iy_per,iz_per, + & ix_rel,iy_rel,iz_rel,ix_dim,iy_dim,iz_dim,x_fact, + & interior,iboundary) + + if ((interior.eq.1).and.(iboundary.eq.1)) then + ibase = ibase + 1 + goto 30 + endif + + layer_max = wn_level(iref) - ibase + 1 + + if (layer_min.eq.wn_level(iref)) then + indep_field = 1 + else + indep_field = 0 + endif + else + layer_max = wn_level(iref) + indep_field = 1 + endif + + if (idebug.eq.1) then + print*,'Return 3' + print*,'layer_min,layer_max,indep_field',layer_min,layer_max,indep_field + print*,'interior,iboundary',interior,iboundary + print*,'ibase = ',ibase + print*,'iref,nref,wn_level(iref)',iref,nref,wn_level(iref) + endif + + + return + + end + + + + +c The function takes a given basis function specified by a corner ixc,iyc,izc +c and a size isz at level wn_c in the oct-tree and returns two integer values. +c (i) interior: +c Value 1 if the basis function is completely within the given +c refinement. +c +c Value 0 if the basis function is without the refinement, or +c overlaps the edges of the refinement, or the edges of the +c primary white noise patch. +c +c (ii) iboundary: +c Value 1 if the basis function is sufficiently far from the +c refinement boundary. +c +c Value 0 otherwise. +c The given refinement is defined at level wn_r in the oct-tree and by the variables +c (ix_rel,iy_rel,iz_rel) which give the location of the refinement relative to +c corner of the white noise patch, (ix_per,iy_per,iz_per) which define the +c periodicity of the white noise patch, and (ix_dim,iy_dim,iz_dim) which +c define the size of the refinement. +c +c +c +c================================================================================= + recursive subroutine inref(ixc,iyc,izc,isz,ir1,ir2,nref,wn_level, + & ix_abs,iy_abs,iz_abs,ix_per,iy_per,iz_per, + & ix_rel,iy_rel,iz_rel,ix_dim,iy_dim,iz_dim,x_fact, + & interior,iboundary) + implicit none + + integer nref + integer ixc,iyc,izc,isz,ir1,ir2 + integer wn_level(nref) + integer*8 ix_abs(nref),iy_abs(nref),iz_abs(nref) + integer*8 ix_per(nref),iy_per(nref),iz_per(nref) + integer*8 ix_rel(nref),iy_rel(nref),iz_rel(nref) + integer*8 ix_dim(nref),iy_dim(nref),iz_dim(nref) + integer interior, iboundary + integer x_fact + + integer*8 ixco,iyco,izco,isize + integer*8 ixref0,iyref0,izref0 + integer*8 ixref1,iyref1,izref1 + integer*8 idist + + integer delta_wn + +c Error checking + if (ir2.lt.ir1) stop 'ir2 - 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 (rcoarse); - - fftw_real *rfine = new fftw_real[ nxf * nyf * nzfp]; - fftw_complex *cfine = reinterpret_cast (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 +#include +#include + +#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("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("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("setup", "shift_x"); + coordinate_system_shift_[1] = -pcf_->getValue("setup", "shift_y"); + coordinate_system_shift_[2] = -pcf_->getValue("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 &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(field); + fftw_complex *cfield = reinterpret_cast(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(field); + fftw_complex *cfield = reinterpret_cast(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 +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 &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(pr0); + pc1 = reinterpret_cast(pr1); + pc2 = reinterpret_cast(pr2); + pc3 = reinterpret_cast(pr3); + pc4 = reinterpret_cast(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< %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< = %g | var(p) = %g", level, sum, sum2); +} + +namespace { +RNG_plugin_creator_concrete creator("PANPHASIA"); +} + +#endif