diff --git a/CMakeLists.txt b/CMakeLists.txt index a82b136..b022c05 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,44 +1,55 @@ -cmake_minimum_required(VERSION 3.10) +cmake_minimum_required(VERSION 2.8.8) project(FENNECS Fortran C) +if(POLICY CMP0074) + cmake_policy(SET CMP0074 NEW) +endif() + find_package(MPI COMPONENTS Fortran REQUIRED) find_package(OpenMP COMPONENTS Fortran C REQUIRED) -find_package(HDF5 COMPONENTS Fortran REQUIRED) +find_package(HDF5 COMPONENTS Fortran C REQUIRED) +#mark_as_advanced_prefix(HDF5_Fortran) +#find_package(MPI REQUIRED) + set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake") +include(CMakeFlagsHandling) +# Compiler flags for debug/optimization +if(${CMAKE_Fortran_COMPILER_ID} MATCHES "Intel") + set(CMAKE_AR ${XIAR}) + add_flags(LANG Fortran TYPE DEBUG -traceback "-check bounds" "-warn unused" -fpp -qopenmp -fpic) + add_flags(LANG Fortran TYPE RELEASE -xHost -O3 -g -warn -fpp -qopenmp -fpic) +elseif(${CMAKE_Fortran_COMPILER_ID} MATCHES "GNU") + add_flags(LANG Fortran TYPE DEBUG -fbounds-check -g -fbacktrace -fpic) + add_flags(LANG Fortran "-fopenmp" "-cpp" "-ffree-line-length-none" "-fpic" "-std=legacy" "-fallow-argument-mismatch") +endif() + # Search and load the FUTILS configuration file if(NOT TARGET futils) - find_package(futils PATHS ${FUTILS}/lib/cmake REQUIRED) + find_package(futils PATHS ${futils_DIR}/lib/cmake REQUIRED) endif() # Search and load the bsplines configuration file if(NOT TARGET bsplines) - find_package(bsplines PATHS ${BSPLINES}/lib/cmake REQUIRED) + find_package(bsplines PATHS ${bsplines_DIR}/lib/cmake REQUIRED) endif() #Search and load the bsplines configuration file if(NOT TARGET sisl) - find_package(sisl REQUIRED) + find_package(sisl MODULE REQUIRED) endif() if(NOT TARGET forsSISL) - find_package(forSISL REQUIRED) + find_package(forSISL MODULE REQUIRED) endif() -include(CMakeFlagsHandling) -# Compiler flags for debug/optimization -if(${CMAKE_Fortran_COMPILER_ID} MATCHES "Intel") - set(CMAKE_AR ${XIAR}) - add_flags(LANG Fortran TYPE DEBUG -traceback "-check bounds" "-warn unused") - add_flags(LANG Fortran TYPE RELEASE -xHost -O3 -g -warn) -elseif(${CMAKE_Fortran_COMPILER_ID} MATCHES "GNU") - add_flags(LANG Fortran TYPE DEBUG -fbounds-check -fbacktrace) +if(POLICY CMP0028) + cmake_policy(SET CMP0028 OLD) endif() - add_subdirectory(src) diff --git a/cmake/Findbsplines.cmake b/cmake/Findbsplines.cmake deleted file mode 100644 index e5b71a7..0000000 --- a/cmake/Findbsplines.cmake +++ /dev/null @@ -1,8 +0,0 @@ -find_library(bsplines_LIBRARY NAMES bsplines fft pppack pputils2 PATHS ${BSPLINES}/lib) -find_path(bsplines_INCLUDES bsplines.mod ${BSPLINES}/include) -find_path(bsplines_INCLUDE_DIR bsplines.mod ${BSPLINES}/include) - -mark_as_advanced(bsplines_LIBRARY bsplines_INCLUDES) - -include(FindPackageHandleStandardArgs) -find_package_handle_standard_args(bsplines DEFAULT_MSG bsplines_LIBRARY) diff --git a/dependencies/random/random_mod.f90 b/dependencies/random/random_mod.f90 index 0baa71e..d6ab698 100644 --- a/dependencies/random/random_mod.f90 +++ b/dependencies/random/random_mod.f90 @@ -1,21 +1,506 @@ -MODULE random +MODULE random_mod ! ! Random module for the random number generator (RNG) ! ! This module, together with the file "random.h", allow to use ! a high-quality, fast, parallel RNG developed by Charles Karney in Princeton University INCLUDE "random.h" CHARACTER(ran_c) :: random_seed_str = '2.71828180' INTEGER, allocatable, DIMENSION(:,:) :: seed ! nbprocs*ran_s size INTEGER, PARAMETER :: n_rand = 1 INTEGER, allocatable:: ran_index(:) ! nbporcs size DOUBLE PRECISION, allocatable :: ran_array(:,:) ! nbprocs*ran_k size -END MODULE random + CONTAINS + + function random(p,r) + implicit none + DOUBLE PRECISION ulp2 + parameter(ulp2=2.0d0**(-47-1)) + REAL ulps + DOUBLE PRECISION mult + parameter(ulps=2.0**(-23),mult=2.0d0**23) + DOUBLE PRECISION random + REAL srandom + integer p + DOUBLE PRECISION r(0:100-1) + if(p.GE.100)then + call rand_batch(p,r(0)) + end if + random=r(p)+ulp2 + p=p+1 + return + entry srandom(p,r) + if(p.GE.100)then + call rand_batch(p,r(0)) + end if + srandom=(int(mult*r(p))+0.5)*ulps + p=p+1 + return + end + + subroutine random_array(y,n,p,r) + implicit none + DOUBLE PRECISION ulp2 + parameter(ulp2=2.0d0**(-47-1)) + REAL ulps + DOUBLE PRECISION mult + parameter(ulps=2.0**(-23),mult=2.0d0**23) + integer n + DOUBLE PRECISION y(0:n-1) + REAL ys(0:n-1) + integer p + DOUBLE PRECISION r(0:100-1) + integer i,k,j + if(n.LE.0)return + k=min(n,100-p) + !$OMP SIMD + do i=0,k-1 + y(i)=r(i+p)+ulp2 + end do + !$OMP END SIMD + p=p+(k) + do j=k,n-1,100 + call rand_batch(p,r(0)) + !$OMP SIMD + do i=j,min(j+100,n)-1 + y(i)=r(i-j+p)+ulp2 + end do + !$OMP END SIMD + p=p+(min(100,n-j)) + end do + return + entry srandom_array(ys,n,p,r) + if(n.LE.0)return + k=min(n,100-p) + do i=0,k-1 + ys(i)=(int(mult*r(i+p))+0.5)*ulps + end do + p=p+(k) + do j=k,n-1,100 + call rand_batch(p,r(0)) + do i=j,min(j+100,n)-1 + ys(i)=(int(mult*r(i-j+p))+0.5)*ulps + end do + p=p+(min(100,n-j)) + end do + return + end + + subroutine rand_batch(p,r) + implicit none + integer p + DOUBLE PRECISION r(0:100-1) + integer i + DOUBLE PRECISION w(0:1009-100-1) + DOUBLE PRECISION tmp + do i=0,63-1 + tmp=r(i)+r(i+100-63) + w(i)=tmp-int(tmp) + end do + do i=63,100-1 + tmp=r(i)+w(i-63) + w(i)=tmp-int(tmp) + end do + do i=100,1009-100-1 + tmp=w(i-100)+w(i-63) + w(i)=tmp-int(tmp) + end do + do i=1009-100,1009-100+63-1 + tmp=w(i-100)+w(i-63) + r(i-1009+100)=tmp-int(tmp) + end do + do i=1009-100+63,1009-1 + tmp=w(i-100)+r(i-1009+100-63) + r(i-1009+100)=tmp-int(tmp) + end do + p=0 + return + end + + subroutine random_init(seed,p,r) + implicit none + integer b + DOUBLE PRECISION del,ulp + parameter(del=2.0d0**(-14),ulp=2.0d0**(-47),b=2**14) + integer seed(0:8-1) + integer p + DOUBLE PRECISION r(0:100-1) + integer i,j,t,s(0:8-1) + logical even + integer z(0:8-1) + integer a0,a1,a2,a3,a4,a5,a6,c0 + data a0,a1,a2,a3,a4,a5,a6,c0/15661,678,724,5245,13656,11852,29,1/ + do i=0,8-1 + s(i)=seed(i) + end do + even=.TRUE. + even=even.AND.(mod(s(7),2).EQ.0) + r(0)=(((s(7)*del+s(6))*del+s(5))*del+int(s(4)/512))*512*del + do j=1,100-1 + t=c0+a0*s(0) + z(0)=mod(t,b) + t=int(t/b)+a0*s(1)+a1*s(0) + z(1)=mod(t,b) + t=int(t/b)+a0*s(2)+a1*s(1)+a2*s(0) + z(2)=mod(t,b) + t=int(t/b)+a0*s(3)+a1*s(2)+a2*s(1)+a3*s(0) + z(3)=mod(t,b) + t=int(t/b)+a0*s(4)+a1*s(3)+a2*s(2)+a3*s(1)+a4*s(0) + z(4)=mod(t,b) + t=int(t/b)+a0*s(5)+a1*s(4)+a2*s(3)+a3*s(2)+a4*s(1)+a5*s(0) + z(5)=mod(t,b) + t=int(t/b)+a0*s(6)+a1*s(5)+a2*s(4)+a3*s(3)+a4*s(2)+a5*s(1)+a6*s(0) + z(6)=mod(t,b) + t=int(t/b)+a0*s(7)+a1*s(6)+a2*s(5)+a3*s(4)+a4*s(3)+a5*s(2)+a6*s(1) + z(7)=mod(t,b) + do i=0,8-1 + s(i)=z(i) + end do + even=even.AND.(mod(s(7),2).EQ.0) + r(j)=(((s(7)*del+s(6))*del+s(5))*del+int(s(4)/512))*512*del + end do + if(even)then + t=c0+a0*s(0) + z(0)=mod(t,b) + t=int(t/b)+a0*s(1)+a1*s(0) + z(1)=mod(t,b) + t=int(t/b)+a0*s(2)+a1*s(1)+a2*s(0) + z(2)=mod(t,b) + t=int(t/b)+a0*s(3)+a1*s(2)+a2*s(1)+a3*s(0) + z(3)=mod(t,b) + t=int(t/b)+a0*s(4)+a1*s(3)+a2*s(2)+a3*s(1)+a4*s(0) + z(4)=mod(t,b) + t=int(t/b)+a0*s(5)+a1*s(4)+a2*s(3)+a3*s(2)+a4*s(1)+a5*s(0) + z(5)=mod(t,b) + t=int(t/b)+a0*s(6)+a1*s(5)+a2*s(4)+a3*s(3)+a4*s(2)+a5*s(1)+a6*s(0) + z(6)=mod(t,b) + t=int(t/b)+a0*s(7)+a1*s(6)+a2*s(5)+a3*s(4)+a4*s(3)+a5*s(2)+a6*s(1) + z(7)=mod(t,b) + do i=0,8-1 + s(i)=z(i) + end do + j=int((s(8-1)*100)/b) + r(j)=r(j)+(ulp) + end if + p=100 + return + end + + subroutine decimal_to_seed(decimal,seed) + implicit none + integer b + parameter(b=2**14) + character*(*)decimal + integer seed(0:8-1) + integer i,ten(0:8-1),c(0:8-1),ch + data ten/10,7*0/ + do i=0,8-1 + seed(i)=0 + c(i)=0 + end do + do i=1,len(decimal) + ch=ichar(decimal(i:i)) + if(ch.GE.ichar('0').AND.ch.LE.ichar('9'))then + c(0)=ch-ichar('0') + call rand_axc(ten,seed,c) + end if + end do + return + end + + subroutine string_to_seed(string,seed) + implicit none + integer b + parameter(b=2**14) + character*(*)string + integer seed(0:8-1) + integer t,i,k,unity(0:8-1),c(0:8-1),ch + data unity/1,7*0/ + do i=0,8-1 + seed(i)=0 + c(i)=0 + end do + do i=1,len(string) + ch=ichar(string(i:i)) + if(ch.GT.ichar(' ').AND.ch.LT.127)then + t=mod(seed(0),2)*(b/2) + do k=0,8-1 + seed(k)=int(seed(k)/2) + if(k.LT.8-1)then + seed(k)=seed(k)+(mod(seed(k+1),2)*(b/2)) + else + seed(k)=seed(k)+(t) + end if + end do + c(0)=ch + call rand_axc(unity,seed,c) + end if + end do + return + end + + subroutine set_random_seed(time,seed) + implicit none + integer time(8) + integer seed(0:8-1) + character*21 c + write(c(1:8),'(i4.4,2i2.2)')time(1),time(2),time(3) + write(c(9:12),'(i1.1,i3.3)') (1-sign(1,time(4)))/2,abs(time(4)) + write(c(13:21),'(3i2.2,i3.3)')time(5),time(6),time(7),time(8) + call decimal_to_seed(c,seed) + return + end + + subroutine seed_to_decimal(seed,decimal) + implicit none + integer pow,decbase,b + parameter(pow=4,decbase=10**pow,b=2**14) + character*(*)decimal + integer seed(0:8-1) + integer z(0:8-1),i,t,j,k + character*36 str + k=-1 + do i=0,8-1 + z(i)=seed(i) + if(z(i).GT.0)k=i + end do + str=' ' + i=9 +90000 continue + i=i-1 + t=0 + do j=k,0,-1 + z(j)=z(j)+t*b + t=mod(z(j),decbase) + z(j)=int(z(j)/decbase) + end do + if(z(max(0,k)).EQ.0)k=k-1 + if(k.GE.0)then + write(str(pow*i+1:pow*(i+1)),'(i4.4)')t + else + write(str(pow*i+1:pow*(i+1)),'(i4)')t + end if + if(k.GE.0)goto 90000 + if(len(decimal).GT.len(str))then + decimal(:len(decimal)-len(str))=' ' + decimal(len(decimal)-len(str)+1:)=str + else + decimal=str(len(str)-len(decimal)+1:) + end if + return + end + + subroutine rand_next_seed(n,ax,cx,y) + implicit none + integer n,ax(0:8-1),cx(0:8-1) + integer y(0:8-1) + integer a(0:8-1),c(0:8-1),z(0:8-1),t(0:8-1),m,i + data z/8*0/ + if(n.EQ.0)return + m=n + do i=0,8-1 + a(i)=ax(i) + c(i)=cx(i) + end do +90000 continue + if(mod(m,2).GT.0)then + call rand_axc(a,y,c) + end if + m=int(m/2) + if(m.EQ.0)return + do i=0,8-1 + t(i)=c(i) + end do + call rand_axc(a,c,t) + do i=0,8-1 + t(i)=a(i) + end do + call rand_axc(t,a,z) + goto 90000 + end + + subroutine next_seed3(n0,n1,n2,seed) + implicit none + integer n0,n1,n2 + integer seed(0:8-1) + integer af0(0:8-1),cf0(0:8-1) + integer ab0(0:8-1),cb0(0:8-1) + integer af1(0:8-1),cf1(0:8-1) + integer ab1(0:8-1),cb1(0:8-1) + integer af2(0:8-1),cf2(0:8-1) + integer ab2(0:8-1),cb2(0:8-1) + data af0/15741,8689,9280,4732,12011,7130,6824,12302/ + data cf0/16317,10266,1198,331,10769,8310,2779,13880/ + data ab0/9173,9894,15203,15379,7981,2280,8071,429/ + data cb0/8383,3616,597,12724,15663,9639,187,4866/ + data af1/8405,4808,3603,6718,13766,9243,10375,12108/ + data cf1/13951,7170,9039,11206,8706,14101,1864,15191/ + data ab1/6269,3240,9759,7130,15320,14399,3675,1380/ + data cb1/15357,5843,6205,16275,8838,12132,2198,10330/ + data af2/445,10754,1869,6593,385,12498,14501,7383/ + data cf2/2285,8057,3864,10235,1805,10614,9615,15522/ + data ab2/405,4903,2746,1477,3263,13564,8139,2362/ + data cb2/8463,575,5876,2220,4924,1701,9060,5639/ + if(n2.GT.0)then + call rand_next_seed(n2,af2,cf2,seed) + else if(n2.LT.0)then + call rand_next_seed(-n2,ab2,cb2,seed) + end if + if(n1.GT.0)then + call rand_next_seed(n1,af1,cf1,seed) + else if(n1.LT.0)then + call rand_next_seed(-n1,ab1,cb1,seed) + end if + entry next_seed(n0,seed) + if(n0.GT.0)then + call rand_next_seed(n0,af0,cf0,seed) + else if(n0.LT.0)then + call rand_next_seed(-n0,ab0,cb0,seed) + end if + return + end + + subroutine rand_axc(a,x,c) + implicit none + integer b + parameter(b=2**14) + integer a(0:8-1),c(0:8-1) + integer x(0:8-1) + integer z(0:8-1),i,j,t + t=0 + do i=0,8-1 + t=c(i)+int(t/b) + do j=0,i + t=t+(a(j)*x(i-j)) + end do + z(i)=mod(t,b) + end do + do i=0,8-1 + x(i)=z(i) + end do + return + end + + subroutine random_gauss(y,n,p,r) + implicit none + integer p + DOUBLE PRECISION r(0:100-1) + integer n + DOUBLE PRECISION y(0:n-1) + integer i + DOUBLE PRECISION pi,theta,z + !qDOUBLE PRECISION random + data pi/3.14159265358979323846264338328d0/ + REAL ys(0:n-1) + REAL spi,stheta,sz + !REAL srandom + data spi/3.14159265358979323846264338328/ + if(n.LE.0)return + call random_array(y,n,p,r(0)) + do i=0,int(n/2)*2-1,2 + theta=pi*(2.0d0*y(i)-1.0d0) + z=sqrt(-2.0d0*log(y(i+1))) + y(i)=z*cos(theta) + y(i+1)=z*sin(theta) + end do + if(mod(n,2).EQ.0)return + theta=pi*(2.0d0*y(n-1)-1.0d0) + z=sqrt(-2.0d0*log(random(p,r(0)))) + y(n-1)=z*cos(theta) + return + entry srandom_gauss(ys,n,p,r) + if(n.LE.0)return + call srandom_array(ys,n,p,r(0)) + do i=0,int(n/2)*2-1,2 + stheta=spi*(2.0*ys(i)-1.0) + sz=sqrt(-2.0*log(ys(i+1))) + ys(i)=sz*cos(stheta) + ys(i+1)=sz*sin(stheta) + end do + if(mod(n,2).EQ.0)return + stheta=spi*(2.0*ys(n-1)-1.0) + sz=sqrt(-2.0*srandom(p,r(0))) + ys(n-1)=sz*cos(stheta) + return + end + + subroutine random_isodist(v,n,p,r) + implicit none + integer p + DOUBLE PRECISION r(0:100-1) + integer n + DOUBLE PRECISION v(0:3*n-1) + integer i + DOUBLE PRECISION pi,costheta,phi + data pi/3.14159265358979323846264338328d0/ + REAL vs(0:3*n-1) + REAL spi,scostheta,sphi + data spi/3.14159265358979323846264338328/ + if(n.LE.0)return + call random_array(v(n),2*n,p,r(0)) + do i=0,n-1 + costheta=2.0d0*v(n+2*i)-1.0d0 + phi=pi*(2.0d0*v(n+2*i+1)-1.0d0) + v(3*i)=cos(phi)*sqrt(1.0d0-costheta**2) + v(3*i+1)=sin(phi)*sqrt(1.0d0-costheta**2) + v(3*i+2)=costheta + end do + return + entry srandom_isodist(vs,n,p,r) + if(n.LE.0)return + call srandom_array(vs(n),2*n,p,r(0)) + do i=0,n-1 + scostheta=2.0*vs(n+2*i)-1.0 + sphi=spi*(2.0*vs(n+2*i+1)-1.0) + vs(3*i)=cos(sphi)*sqrt(1.0-scostheta**2) + vs(3*i+1)=sin(sphi)*sqrt(1.0-scostheta**2) + vs(3*i+2)=scostheta + end do + return + end + + subroutine random_cosdist(v,n,p,r) + implicit none + integer p + DOUBLE PRECISION r(0:100-1) + integer n + DOUBLE PRECISION v(0:3*n-1) + integer i + DOUBLE PRECISION pi,costheta2,phi + data pi/3.14159265358979323846264338328d0/ + REAL vs(0:2*n-1) + REAL spi,scostheta2,sphi + data spi/3.14159265358979323846264338328/ + if(n.LE.0)return + call random_array(v(n),2*n,p,r(0)) + do i=0,n-1 + costheta2=v(n+2*i) + phi=pi*(2.0d0*v(n+2*i+1)-1.0d0) + v(3*i)=cos(phi)*sqrt(1.0d0-costheta2) + v(3*i+1)=sin(phi)*sqrt(1.0d0-costheta2) + v(3*i+2)=sqrt(costheta2) + end do + return + entry srandom_cosdist(vs,n,p,r) + if(n.LE.0)return + call srandom_array(vs(n),2*n,p,r(0)) + do i=0,n-1 + scostheta2=vs(n+2*i) + sphi=spi*(2.0*vs(n+2*i+1)-1.0) + vs(3*i)=cos(sphi)*sqrt(1.0-scostheta2) + vs(3*i+1)=sin(sphi)*sqrt(1.0-scostheta2) + vs(3*i+2)=sqrt(scostheta2) + end do + return + end + + + END MODULE random_mod diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 7ab91b5..7b6025e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,83 +1,82 @@ project(fennecs_src Fortran C) add_executable(fennecs) set(SRCS main.f90 basic_mod.f90 newrun.f90 restart.f90 auxval.f90 inital.f90 resume.f90 start.f90 diagnose.f90 stepon.f90 tesend.f90 endrun.f90 chkrst.f90 mv2bk.f90 constants.f90 fields_mod.f90 beam_mod.f90 mpihelper_mod.f90 sort_mod.f90 distrib_mod.f90 maxwsrce_mod.f90 celldiag_mod.f90 geometry_mod.f90 -random_mod.f90 neutcol_mod.f90 particletypes_mod.f90 splinebound_mod.f90 weighttypes_mod.f90 psupply_mod.f90 ion_induced_mod.f90 -incomplete_gamma_mod.f90 materials_mod.f90 extra.c -random.f -randother.f -elliptic_mod.f90 magnet_mod.f90 ../dependencies/random/random_mod.f90 -../dependencies/random/random.f -../dependencies/random/randother.f ../dependencies/incomplete_gamma/incomplete_gamma_mod.f90 ../dependencies/elliptic/elliptic_mod.f90 ) target_sources(fennecs PRIVATE ${SRCS}) -set_property(SOURCE ${SRCS} APPEND PROPERTY COMPILE_OPTIONS -cpp -fpp -qopenmp ) +set_property(SOURCE ${SRCS} APPEND PROPERTY COMPILE_OPTIONS -cpp -fopenmp ) -include_directories(SYSTEM ${MPI_INCLUDE_PATH} ${forSISL_INCLUDE_DIR} ${bsplines_INCLUDES} ${futils_INCLUDES}) +include_directories(SYSTEM ${MPI_INCLUDE_PATH} ${forSISL_INCLUDE_DIR} ${BSPLINES_MODS} ${FUTILS_MODS}) if(MKL_Fortran_FLAGS) separate_arguments(MKL_Fortran_FLAGS) target_compile_options(fennecs PUBLIC ${MKL_Fortran_FLAGS}) target_link_options(fennecs PUBLIC ${MKL_Fortran_FLAGS}) endif() add_custom_command(TARGET fennecs POST_BUILD COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_BINARY_DIR}/fennecs ${CMAKE_CURRENT_BINARY_DIR}/../fennecs) target_include_directories(fennecs PRIVATE $ - futils) - -target_link_libraries(fennecs PUBLIC futils bsplines MPI::MPI_Fortran -OpenMP::OpenMP_Fortran -${BLAS_LIBRARIES} -${MUMPS_LIBRARIES} -${LAPACK_LIBRARIES} -${forSISL_LIBRARY} -${sisl_LIBRARY} -${futils_LIBRARY} -${bsplines_LIBRARY} -) + ${futils_INCLUDE_DIR}) +target_link_libraries(fennecs PUBLIC MPI::MPI_Fortran + futils bsplines + OpenMP::OpenMP_Fortran + #${BLAS_LIBRARIES} + #${MUMPS_LIBRARIES} + #${LAPACK_LIBRARIES} + ${forSISL_LIBRARY} + ${sisl_LIBRARY} + ${HDF5_Fortran_LIBRARIES} + ) +get_cmake_property(_variableNames VARIABLES) +list (SORT _variableNames) +foreach (_variableName ${_variableNames}) + message(STATUS "${_variableName}=${${_variableName}}") +endforeach() +get_target_property(OUT fennecs LINK_LIBRARIES) +message(STATUS ${OUT}) diff --git a/src/auxval.f90 b/src/auxval.f90 index 7115220..f5634df 100644 --- a/src/auxval.f90 +++ b/src/auxval.f90 @@ -1,46 +1,46 @@ SUBROUTINE auxval ! USE constants USE basic USE fields USE beam - USE random + USE random_mod USE omp_lib ! ! Set auxiliary values ! IMPLICIT NONE ! !________________________________________________________________________________ IF(mpirank .eq. 0) WRITE(*,'(a/)') '=== Set auxiliary values ===' !________________________________________________________________________________ ! ! Define Zindex bounds for the MPI process in the case of sequential run ALLOCATE( Zbounds(0:mpisize), Nplocs_all(0:mpisize-1)) Zbounds=0 Zbounds(mpisize) = nz IF(mpisize.gt.1) THEN ! If we use mpi we define the left and right nodes used for communications leftproc=mpirank-1 IF(mpirank .eq. 0 .AND. partperiodic) THEN leftproc = mpisize-1 ELSE IF (mpirank .eq. 0) THEN leftproc=-1 END IF rightproc = mpirank+1 IF(mpirank .eq. mpisize-1 .AND. partperiodic) THEN rightproc = 0 ELSE IF (mpirank .eq. mpisize-1) THEN rightproc=-1 END IF ELSE leftproc=-1 rightproc=-1 END IF WRITE(*,*) "mpirank, left, right", mpirank, leftproc, rightproc !________________________________________________________________________________ !________________________________________________________________________________ END SUBROUTINE auxval diff --git a/src/basic_mod.f90 b/src/basic_mod.f90 index 98654f3..fc01a77 100644 --- a/src/basic_mod.f90 +++ b/src/basic_mod.f90 @@ -1,450 +1,450 @@ MODULE basic ! USE hashtable USE constants USE bsplines USE mumps_bsplines USE futils USE mpihelper - use random + use random_mod IMPLICIT NONE ! ! Basic module for time dependent problems ! CHARACTER(len=128) :: label1, label2, label3, label4 ! ! BASIC Namelist ! LOGICAL :: nlres = .FALSE. !< Restart flag LOGICAL :: nlsave = .TRUE. !< Checkpoint (save) flag LOGICAL :: newres=.FALSE. !< New result HDF5 file LOGICAL :: nlxg=.FALSE. !< Show graphical interface Xgrafix LOGICAL :: nlmaxwellsource = .FALSE. !< Activate the maxwell source INTEGER :: nrun=1 !< Number of time steps to run REAL(kind=db) :: job_time=3600.0 !< Time allocated to this job in seconds REAL(kind=db) :: tmax=100000.0 !< Maximum simulation time REAL(kind=db) :: extra_time=60.0 !< Extra time allocated REAL(kind=db) :: dt=1 !< Time step REAL(kind=db) :: time=0 !< Current simulation time (Init from restart file) ! ! Other basic global vars and arrays ! INTEGER :: jobnum !< Job number INTEGER :: step=0 !< Calculation step of this run INTEGER :: cstep=0 !< Current step number (Init from restart file) LOGICAL :: nlend=.false. !< Signal end of run INTEGER :: ierr !< Integer used for MPI INTEGER :: it0d=1 !< Number of iterations between 0d values writes to hdf5 INTEGER :: it2d=100 !< Number of iterations between 2d values writes to hdf5 INTEGER :: itparts=1000 !< Number of iterations between particles values writes to hdf5 INTEGER :: ittext=10 !< Number of iterations between text outputs in the console INTEGER :: itrestart=10000 !< Number of iterations between save of restart.h5 file INTEGER :: ittracer=100 !< Number of iterations between save of traced particles position and velocity INTEGER :: itcelldiag=100000 !< Number of iterations between save of celldiag diagnostic INTEGER :: nbcelldiag=0 !< Number of celldiagnostics INTEGER :: itgraph !< Number of iterations between graphical interface updates INTEGER :: mpirank !< MPIrank of the current processus INTEGER :: mpisize !< Size of the MPI_COMM_WORLD communicator INTEGER :: rightproc !< Rank of next processor in the z decomposition INTEGER :: leftproc !< Rank of previous processor in the z decomposition ! ! List of logical file units INTEGER :: lu_in = 90 !< File duplicated from STDIN INTEGER :: lu_stop = 91 !< stop file, see subroutine TESEND INTEGER :: lu_partfile = 120 !< particle loading file, see beam::loadpartfile ! ! HDF5 file CHARACTER(len=256) :: resfile = "results.h5" !< Main result file CHARACTER(len=256) :: rstfile = "restart.h5" !< Restart file CHARACTER(len=256) :: magnetfile = "" !< H5 file containing the magnetic field definition where r,z are in m and Br, Bz are in T CHARACTER(len=256) :: partfile(10)="" !< Particle loading file CHARACTER(len=256) :: addedtestspecfile(10)="" !< Particle file list for added particles at restart INTEGER :: fidres !< File ID for resfile INTEGER :: fidrst !< File ID for restart file TYPE(BUFFER_TYPE) :: hbuf0 !< Hashtable for 0d var ! ! Plasma parameters LOGICAL :: nlclassical= .FALSE. !< If true, solves the equation of motion according to classical !! dynamics LOGICAL :: partperiodic= .TRUE. !< Sets if the particles boundary conditions are periodic or open INTEGER :: nbaddtestspecies=0 !< On restart number of files to read to add test particles INTEGER :: nplasma !< Number of macro-particles on initialisation REAL(kind=db) :: qsim !< Charge of superparticles [C] REAL(kind=db) :: msim !< Mass of superparticles [kg] REAL(kind=db) :: partmass=me !< Mass of physical particle [kg] INTEGER :: nbspecies = 1 !< Number of particles species also counting tracing particles INTEGER :: npartsalloc = 0 !< Size of particle memory allocated at the begining of the simulation INTEGER :: nblock !< Number of slices in Z for stable distribution initialisation REAL(kind=db) :: n0 !< Physical plasma density parameter [m-3] used in distribtype=1 and for time scales normlisation REAL(kind=db) :: plasmadim(4) !< Zmin Zmax Rmin Rmax values for plasma particle loading [m] REAL(kind=db) :: H0=0 !< Initial value of Hamiltonian for distribution 2 [J] REAL(kind=db) :: P0=0 !< Initial canonical angular momentum for distribution 2 [kg m^2/s] REAL(kind=db) :: temp !< Initial temperature of plasma [K] for distribtype=1 REAL(kind=db) :: weights_scale=1.0 !< Scale factor for the particle weights on restart (only for newres=.true.) REAL(kind=db) :: temprescale = -1.0 !< Factor used for temperature rescaling in case of a restart (<0 -> no rescaling) Currently not implemented INTEGER :: samplefactor =-1 !< Factor used for the up-sampling of the particles number ! ! Fields parameters REAL(kind=db) :: B0 !< Max magnitude of magnetic field [T] and normalisation factor for magnetic field ! If magnetfile ='' The magnetic field is one of a magnetic mirror with maximum amplitude on axis of B0 ! and REAL(kind=db) :: Rcurv = 1.0 !< Magnetic field curvature coefficient REAL(kind=db) :: Width = 1.0 !< Distance between two magnetic mirrors REAL(kind=db), allocatable :: Bz(:), Br(:) !< Normalised magnetic field components REAL(kind=db), allocatable :: Athet(:) !< Theta component of the magnetic vector potential Tm INTEGER :: bscaling = -1 !< if >0 rescale the magnetic field read from h5 file before calculating value at grid points, if <0 rescale after interpolation, if = 0 doesn't rescale INTEGER :: Dimensionality = 1 !< Defines the dimensionality of the solution 1: (r,z) 2: (r,theta) 3: (r,theta,z) LOGICAL :: nlperiod(3)=(/.false.,.false.,.true./)!< Set periodic splines on or off (z,r,theta) LOGICAL :: nlPhis= .TRUE. !< Calculate self consistent electric field flag LOGICAL :: nlfreezephi= .FALSE. !< Freeze the Poisson solver to the field obtained at (re-)start INTEGER :: femorder(3) = 2 !< FEM order (z,r,theta) INTEGER :: ngauss(3) = 6 !< Number of gauss points for FEM integration (z,r,theta) LOGICAL :: nlppform =.TRUE. !< Defines if spline evaluation is done using ppform (faster with true) INTEGER, SAVE :: nrank(3) !< Number of splines in all directions (z,r,theta) TYPE(spline2d), SAVE :: splrz !< Spline at r and z for total electric field TYPE(spline2d), SAVE :: splrz_ext !< Spline at r and z for external electric field TYPE(spline2d), SAVE :: splrthet !< Spline at r and theta for total electric field TYPE(spline2d), SAVE :: splrthet_ext !< Spline at r and theta for external electric field REAL(kind=db) :: potinn=0 !< Electric potential at the inner metallic wall REAL(kind=db) :: potout=0 !< Electric potential at the outer metallic wall REAL(kind=db), allocatable :: Ez(:), Er(:), Ethet(:) !< Normalised electric field components ( ext+self ) REAL(kind=db), allocatable :: pot(:) !< Normalised electrostatic potential ( ext+self ) REAL(kind=db), allocatable :: Ezxt(:), Erxt(:), Ethetxt(:) !< Normalised external Electric field components REAL(kind=db), allocatable :: potxt(:) !< Normalised external Electro static potential REAL(kind=db) :: radii(11) !< Inner and outer radius of cylinder and radii of fine mesh region [m] INTEGER :: distribtype=1 !< Type of distribution function used to load the particles !!1: gaussian, 2: Stable as defined in 4.95 of Davidson, 7 use particle input file REAL(kind=db) :: lz(11) !< Lower and upper cylinder limits in z direction [m] REAL(kind=db), DIMENSION(:), ALLOCATABLE, SAVE:: rhs !< right hand side of the poisson equation solver REAL(kind=db), DIMENSION(:), ALLOCATABLE, SAVE:: volume !< Volume covered by each spline for density calculation INTEGER :: nz !< Total Number of grid intervals in z INTEGER :: nnz(10) !< Number of grid intervals in z INTEGER :: nsubz=10 !< Number of sub-intervals in z INTEGER :: nr !< Total number of grid intervals in r INTEGER :: nnr(10) !< Number of grid intervals in r in each subdomain INTEGER :: nsubr=10 !< Number of sub-intervals in r INTEGER :: nthet=0 !< Total Number of grid intervals in theta REAL(kind=db) :: dz(10) !< Cell size in z REAL(kind=db) :: dr(10) !< Cell size in r for each region Real(kind=db) :: dthet !< Cell size in theta REAL(kind=db) :: invdz(10), invdr(10) !< inverse of the grid cell step Real(kind=db) :: invdthet !< inverse of the theta grid cell step REAL(kind=db), ALLOCATABLE :: zgrid(:) !< Nodes positions in longitudinal direction REAL(kind=db), ALLOCATABLE :: rgrid(:) !< Nodes positions in radial direction REAL(kind=db), ALLOCATABLE :: thetgrid(:) !< Nodes positions in azimuthal direction REAL(kind=db) :: bnorm,enorm,vnorm,tnorm,rnorm,phinorm,qnorm !< Normalization constants REAL(kind=db) :: omegac !< cylotronic frequency at B0 [1/s] REAL(kind=db) :: omegap !< Plasma frequency at n0 [1/s] INTEGER, DIMENSION(:), ALLOCATABLE :: Zbounds !< Index of bounds for local processus in Z direction for MPI decomposition CONTAINS ! !================================================================================ SUBROUTINE basic_data ! ! Define basic data ! use mpihelper USE omp_lib - Use random + Use random_mod IMPLICIT NONE ! ! Local vars and arrays CHARACTER(len=256) :: inputfilename INTEGER :: i, nbprocs ! NAMELIST /BASIC/ job_time, extra_time, nrun, tmax, dt, nlres, nlsave, newres, nlxg, & & nplasma, potinn, potout, B0, lz, n0, nz, nnz, nnr, femorder, ngauss, & & nlppform, plasmadim, radii, temp, Rcurv, width, it0d, it2d, itparts, ittext, & & resfile, rstfile, itgraph, nlPhis, distribtype, nblock, nlclassical, H0, P0, partperiodic, & & temprescale, samplefactor, nlmaxwellsource, npartsalloc, partfile, partmass, nbspecies, & & ittracer, itcelldiag, nbcelldiag, magnetfile, weights_scale, nlfreezephi, nbaddtestspecies, & & addedtestspecfile, bscaling, nthet, Dimensionality !________________________________________________________________________________ ! 1. Process Standard Input File ! IF(COMMAND_ARGUMENT_COUNT().NE.1)THEN WRITE(*,*)'ERROR, ONE COMMAND-LINE ARGUMENT REQUIRED, STOPPING' STOP ENDIF CALL GET_COMMAND_ARGUMENT(1,inputfilename) OPEN(UNIT=lu_in,FILE=trim(inputfilename),ACTION='READ') IF(mpirank .eq. 0) THEN !________________________________________________________________________________ ! 1. Label the run ! READ(lu_in,'(a)') label1 READ(lu_in,'(a)') label2 READ(lu_in,'(a)') label3 READ(lu_in,'(a)') label4 ! WRITE(*,'(12x,a/)') label1(1:len_trim(label1)) WRITE(*,'(12x,a/)') label2(1:len_trim(label2)) WRITE(*,'(12x,a/)') label3(1:len_trim(label3)) WRITE(*,'(12x,a/)') label4(1:len_trim(label4)) !________________________________________________________________________________ ! 2. Read in basic data specific to run ! READ(lu_in,basic) WRITE(*,basic) #if _DEBUG==1 WRITE(*,*) "Compiled in debug mode" #endif ELSE READ(lu_in,basic) END IF CALL mpitypes_init ! initialize all mpi types that will be needed in the simulation WRITE(*,'(a,i4.2,a,i4.2,a)')"Running on ",mpisize," tasks with", omp_get_max_threads() ," openMP threads" WRITE(*,*)"db given kind", db IF(samplefactor .gt. 1 .and. .not. newres) THEN IF(mpirank.eq.0) WRITE(*,*)"To increase the number of particles, you need to create a new result file (set newres to 1)" CALL MPI_abort(MPI_COMM_WORLD,-1,ierr) END IF IF (npartsalloc .lt. nplasma) THEN npartsalloc=nplasma END IF ! Total number of intervals nr=sum(nnr) if (any(nnz.gt.0)) then nz=sum(nnz) else nnz(1)=nz end if ! Normalisation constants if(nplasma .gt. 0) then qsim=pi*(plasmadim(2)-plasmadim(1))*(plasmadim(4)**2-plasmadim(3)**2)*n0*elchar/nplasma else qsim=sign(n0,elchar) end if msim=abs(qsim)/elchar*partmass vnorm=vlight omegac=sign(elchar,qsim)/partmass*B0 omegap=sqrt(elchar**2*abs(n0)/partmass/eps_0) tnorm=min(abs(1/omegac),abs(1/omegap)) rnorm=vnorm*tnorm bnorm=B0 enorm=vlight*bnorm phinorm=enorm*rnorm ! Normalised boundary conditions potinn=potinn/phinorm potout=potout/phinorm ! Normalised dt dt=dt/tnorm ! Characteristic frequencies and normalised volume IF(mpirank .eq. 0) THEN IF(abs(omegap).GT. abs(omegac)) THEN WRITE(*,'(a,3(1pe12.3))') 'omegap, omegac, omegap/omegac', omegap, omegac, omegap/omegac ELSE WRITE(*,'(a,3(1pe12.3))') 'omegap, omegac, omegac/omegap', omegap, omegac, omegac/omegap END IF END IF ! Construction of the mesh rgrid in r and zgrid in z and its normalisation CALL mesh rgrid=rgrid/rnorm zgrid=zgrid/rnorm dz=dz/rnorm dr=dr/rnorm Where(dr.gt.0) invdr=1/dr Where(dz.gt.0) invdz=1/dz !invdz=1/dz ! Initialize random number generator nbprocs = omp_get_max_threads() allocate(seed(ran_s,nbprocs), ran_index(nbprocs), ran_array(ran_k,nbprocs)) IF(.false.) then call date_and_time(time=random_seed_str) CALL MPI_BCAST(random_seed_str,10,MPI_CHARACTER,0,MPI_COMM_WORLD,ierr) write(*,*) "MPI seed:", mpirank, random_seed_str end if Do i=1,nbprocs ! Generate seed from the default seed-string in random module CALL decimal_to_seed(random_seed_str, seed(:,i)) ! Generate a different seed for each processor from the mother seed CALL next_seed(mpirank*nbprocs+i,seed(:,i)) ! Initialize the random array (first hundred numbers) CALL random_init(seed(:,i), ran_index(i), ran_array(:,i)) end do ! END SUBROUTINE basic_data !================================================================================ SUBROUTINE daytim(str) ! ! Print date and time ! IMPLICIT NONE ! CHARACTER(len=*), INTENT(in) :: str ! ! Local vars and arrays CHARACTER(len=16) :: d, t, dat, functime !________________________________________________________________________________ ! CALL DATE_AND_TIME(d,t) dat=d(7:8) // '/' // d(5:6) // '/' // d(1:4) functime=t(1:2) // ':' // t(3:4) // ':' // t(5:10) WRITE(*,'(a,1x,a,1x,a)') str, dat(1:10), functime(1:12) ! END SUBROUTINE daytim !================================================================================ SUBROUTINE timera(cntrl, str, eltime) ! ! Timers (cntrl=0/1 to Init/Update) ! IMPLICIT NONE INTEGER, INTENT(in) :: cntrl CHARACTER(len=*), INTENT(in) :: str REAL(kind=db), OPTIONAL, INTENT(out) :: eltime ! INTEGER, PARAMETER :: ncmax=128 INTEGER, SAVE :: icall=0, nc=0 REAL(kind=db), SAVE :: startt0=0.0 CHARACTER(len=16), SAVE :: which(ncmax) INTEGER :: lstr, found, i REAL(kind=db) :: seconds REAL(kind=db), DIMENSION(ncmax), SAVE :: startt = 0.0, endt = 0.0 !________________________________________________________________________________ IF( icall .EQ. 0 ) THEN icall = icall+1 startt0 = seconds() END IF lstr = LEN_TRIM(str) IF( lstr .GT. 0 ) found = loc(str) !________________________________________________________________________________ ! SELECT CASE (cntrl) ! CASE(-1) ! Current wall time IF( PRESENT(eltime) ) THEN eltime = seconds() - startt0 ELSE WRITE(*,'(/a,a,1pe10.3/)') "++ ", ' Wall time used so far = ', seconds() - startt0 END IF ! CASE(0) ! Init Timer IF( found .EQ. 0 ) THEN ! Called for the 1st time for 'str' nc = nc+1 which(nc) = str(1:lstr) found = nc END IF startt(found) = seconds() ! CASE(1) ! Update timer endt(found) = seconds() - startt(found) IF( PRESENT(eltime) ) THEN eltime = endt(found) ELSE WRITE(*,'(/a,a,1pe10.3/)') "++ "//str, ' wall clock time = ', endt(found) END IF ! CASE(2) ! Update and reset timer endt(found) = endt(found) + seconds() - startt(found) startt(found) = seconds() IF( PRESENT(eltime) ) THEN eltime = endt(found) END IF ! CASE(9) ! Display all timers IF( nc .GT. 0 ) THEN WRITE(*,'(a)') "Timer Summary" WRITE(*,'(a)') "=============" DO i=1,nc WRITE(*,'(a20,2x,2(1pe12.3))') TRIM(which(i))//":", endt(i) END DO END IF ! END SELECT ! CONTAINS INTEGER FUNCTION loc(funcstr) CHARACTER(len=*), INTENT(in) :: funcstr INTEGER :: j, ind loc = 0 DO j=1,nc ind = INDEX(which(j), funcstr(1:lstr)) IF( ind .GT. 0 .AND. LEN_TRIM(which(j)) .EQ. lstr ) THEN loc = j EXIT END IF END DO END FUNCTION loc END SUBROUTINE timera !================================================================================ !--------------------------------------------------------------------------- !> @author !> Patryk Kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Creates the mesh in r, theta and z direction for calculating the electric and magnetic fields. !--------------------------------------------------------------------------- SUBROUTINE mesh INTEGER :: j,i,k ALLOCATE(zgrid(0:nz),rgrid(0:nr),thetgrid(0:nthet)) !dz=(lz(2)-lz(1))/nz k=0 nsubz=count(nnz.gt.0) zgrid(0)=lz(1) do i=1,nsubz dz(i)=(lz(i+1)-lz(i))/nnz(i) if (nnz(i).gt.0) then DO j=1,nnz(i) zgrid(j+k)=lz(i)+j*dz(i) END DO end if k=k+nnz(i) end do nsubr=count(nnr.gt.0) k=0 rgrid(0)=radii(1) do i=1,nsubr dr(i)=(radii(i+1)-radii(i))/nnr(i) if (nnr(i).gt.0) then DO j=1,nnr(i) rgrid(j+k)=radii(i)+j*dr(i) END DO end if k=k+nnr(i) end do thetgrid(0)=0 dthet=0.0_db invdthet=0.0_db if(nthet.gt.0) then dthet=2*pi/nthet invdthet=1/dthet do i=1,nthet thetgrid(i)=i*dthet end do end if END SUBROUTINE mesh END MODULE basic diff --git a/src/distrib_mod.f90 b/src/distrib_mod.f90 index 2993ddf..daf7a22 100644 --- a/src/distrib_mod.f90 +++ b/src/distrib_mod.f90 @@ -1,315 +1,315 @@ !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: distrib ! !> @author !> Guillaume Le Bars EPFL/SPC !> Trach Minh Tran EPFL/SPC ! ! DESCRIPTION: !> Contains the routines used to load several type of distribution functions !------------------------------------------------------------------------------ MODULE distrib USE constants USE omp_lib - USE random + USE random_mod ! IMPLICIT NONE contains !--------------------------------------------------------------------------- !> @author !> Trach Minh Tran EPFL/SPC ! ! DESCRIPTION: !> !> @brief Load a uniform distribution using the Hammersley's sequence (0<=y<=1). !> (nbase = 0 => Random sampling) ! !> @param [in] nbase Base of the Hammersley's sequence. (0 => Random) !> @param [out] y array containing the random variables. !--------------------------------------------------------------------------- SUBROUTINE loduni(nbase, y) ! ! Load a uniform distribution using the Hammersley's sequence. ! (nbase = 0 => Random sampling) ! INTEGER, INTENT(in) :: nbase REAL(kind=db), INTENT(inout) :: y(:) INTEGER :: n, i, ompthread ! n = SIZE(y) ! SELECT CASE (nbase) CASE(0) ompthread=omp_get_thread_num()+1 CALL random_array(y,n,ran_index(ompthread),ran_array(:,ompthread)) !CALL RANDOM_NUMBER(y) CASE(1) DO i=1,n y(i) = (i-0.5_db)/n END DO CASE(2:) DO i=1,n y(i) = rev(nbase,i) END DO CASE default WRITE(*,'(a,i5)') 'Invalid value of NBASE =', nbase END SELECT ! END SUBROUTINE loduni !--------------------------------------------------------------------------- !> @author !> Trach Minh Tran EPFL/SPC ! ! DESCRIPTION: !> !> @brief Load a gaussian distribution using a uniform distribution according to the Hammersley's sequence. !> (nbase = 0 => Random sampling) ! !> @param [in] nbase Base of the Hammersley's sequence. (0 => Random) !> @param [out] y array containing the random variables. !--------------------------------------------------------------------------- SUBROUTINE lodgaus(nbase, y) ! ! Sample y from the Gauss distributrion ! REAL(kind=db), INTENT(inout) :: y(:) INTEGER, INTENT(in) :: nbase ! INTEGER :: n, i, iflag REAL(kind=db) :: r(SIZE(y)) REAL(kind=db) :: t ! REAL(kind=db) :: c0, c1, c2 REAL(kind=db) :: d1, d2, d3 DATA c0, c1, c2/2.515517_db, 0.802853_db, 0.010328_db/ DATA d1, d2, d3/1.432788_db, 0.189269_db, 0.001308_db/ ! n = SIZE(y) CALL loduni(nbase, r) r = 1.0E-5 + 0.99998*r ! DO i=1,n iflag = 1 IF (r(i) .GT. 0.5) THEN r(i) = 1.0 - r(i) iflag = -1 END IF t = SQRT(LOG(1.0/(r(i)*r(i)))) y(i) = t - (c0+c1*t+c2*t**2) / (1.0+d1*t+d2*t**2+d3*t**3) y(i) = y(i) * iflag END DO y = y - SUM(y)/REAL(n) END SUBROUTINE lodgaus !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Load the array y according to a probability distribution function proportional to 1/r !> between ra and rb ! !> @param [in] nbase Base of the Hammersley's sequence. (0 => Random) !> @param [in] ra Lower limit for the values in y !> @param [in] rb Upper limit for the values in y !> @param [out] y array containing the random variables. !--------------------------------------------------------------------------- SUBROUTINE lodinvr(nbase, y, ra, rb) ! ! REAL(kind=db), INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase REAL(kind=db), INTENT(in) :: rb, ra ! REAL(kind=db) :: r(SIZE(y)) ! CALL loduni(nbase, r) r=r*log(rb/ra) y=ra*exp(r) END SUBROUTINE lodinvr !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Load the array y according to a probability distribution function proportional to r !> between ra and rb ! !> @param [in] nbase Base of the Hammersley's sequence. (0 => Random) !> @param [in] ra Lower limit for the values in y !> @param [in] rb Upper limit for the values in y !> @param [out] y array containing the random variables. !--------------------------------------------------------------------------- SUBROUTINE lodlinr(nbase, y, ra, rb) REAL(kind=db), INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase REAL(kind=db), INTENT(in) :: rb, ra ! REAL(kind=db) :: r(SIZE(y)) ! CALL loduni(nbase, r) y=sqrt(r*(rb**2-ra**2)+ra**2) END SUBROUTINE lodlinr !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Load the array y according to a uniform probability distribution function !> between ra and rb ! !> @param [in] nbase Base of the Hammersley's sequence. (0 => Random) !> @param [in] ra Lower limit for the values in y !> @param [in] rb Upper limit for the values in y !> @param [out] y array containing the random variables. !--------------------------------------------------------------------------- SUBROUTINE lodunir(nbase, y, ra, rb) REAL(kind=db), INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase REAL(kind=db), INTENT(in) :: rb, ra ! REAL(kind=db) :: r(SIZE(y)) ! CALL loduni(nbase, r) y=ra+(rb-ra)*r END SUBROUTINE lodunir !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Load the array y according to a gaussian probability distribution function !> around the mean of rb and ra and a sigma of (rb-ra)/10 ! !> @param [in] nbase Base of the Hammersley's sequence. (0 => Random) !> @param [in] ra Lower limit for the values in y !> @param [in] rb Upper limit for the values in y !> @param [out] y array containing the random variables. !--------------------------------------------------------------------------- SUBROUTINE lodgausr(nbase, y, ra, rb) REAL(kind=db), INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase REAL(kind=db), INTENT(in) :: rb, ra ! REAL(kind=db) :: r(SIZE(y)) ! CALL lodgaus(nbase, r) y=0.5*(rb+ra)+ 0.1*(rb-ra)*r END SUBROUTINE lodgausr !--------------------------------------------------------------------------- !> @author !> Trach Minh Tran EPFL/SPC ! ! DESCRIPTION: !> !> @brief Return an element of the Hammersley's sequence ! !> @param [in] nbase Base of the Hammersley's sequence. (0 => Random) !> @param [in] i Index of the sequence. !> @param [out] rev Returned element of the Hammersley's sequence !--------------------------------------------------------------------------- REAL(kind=db) FUNCTION rev(nbase,i) ! INTEGER, INTENT(IN) :: nbase ! Base of the Hammersley's sequence (IN) INTEGER, INTENT(IN) :: i ! Index of the sequence (IN) ! ! Local vars INTEGER :: j1, j2 REAL(kind=db) :: xs, xsi !----------------------------------------------------------------------- xs = 0.0 xsi = 1.0 j2 = i DO xsi = xsi/nbase j1 = j2/nbase xs = xs + (j2-nbase*j1)*xsi j2 = j1 IF( j2.LE.0 ) EXIT END DO rev = xs END FUNCTION rev !--------------------------------------------------------------------------- !> @author !> Trach Minh Tran EPFL/SPC ! ! DESCRIPTION: !> !> @brief Modified Bessel functions of the first kind of the first order 1 ! !--------------------------------------------------------------------------- FUNCTION bessi1(x) REAL(kind=db) :: bessi1,x REAL(kind=db) :: ax REAL(kind=db) p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9 /0.39894228d0, -0.3988024d-1, -0.362018d-2, 0.163801d-2,-0.1031555d-1, & & 0.2282967d-1, -0.2895312d-1, 0.1787654d-1,-0.420059d-2/ if (abs(x).lt.3.75) then y=(x/3.75)**2 bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))) else ax=abs(x) y=3.75/ax bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))) if(x.lt.0.)bessi1=-bessi1 endif return END FUNCTION bessi1 ! find the index for the value in array closest but inferior to val function closest(array, val, siz) real(kind=db):: array(1:), val integer:: closest integer:: siz integer:: i, j, mid ! Corner cases if (val.lt.array(1)) then closest=-1 RETURN end if if (val .ge. array(siz)) then closest=siz return end if !Start binary search i=1 j=siz mid=1 do while (i .lt. j) mid = (i + j) / 2 if (val .ge. array(mid) .and. val .lt. array(mid+1)) then closest = mid return end if ! if val > array(mid) search in right if (val .gt. array(mid)) then ! update i i = mid+1 else ! if val < array(mid) search in left ! update j j = mid end if end do closest=-1 ! not found end function closest END MODULE distrib \ No newline at end of file diff --git a/src/fields_mod.f90 b/src/fields_mod.f90 index ae3ab3f..8706b10 100644 --- a/src/fields_mod.f90 +++ b/src/fields_mod.f90 @@ -1,1559 +1,1558 @@ !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: fields ! !> @author !> Patryk Kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> Module responsible for initializing the magnetic field, solving the Poisson equation and computing the moments of the particles distribution function !------------------------------------------------------------------------------ MODULE fields USE constants USE basic, ONLY: nr, nz, zgrid, rgrid, Br, Bz, Er, Ez, femorder, ngauss, nlppform, pot, Athet, & & splrz, splrz_ext, nlperiod, phinorm, nlPhis, nrank, mpirank, mpisize, step, it2d, timera, potxt, erxt, ezxt, splrthet,splrthet_ext USE beam, ONLY: partslist USE bsplines USE mumps_bsplines use mpi Use omp_lib Use mpihelper, ONLY: db_type USE particletypes IMPLICIT NONE REAL(kind=db), allocatable, SAVE :: matcoef(:, :), phi_spline(:), vec1(:), vec2(:) REAL(kind=db), allocatable, SAVE :: loc_moments(:, :), loc_rhs(:),omp_loc_rhs(:,:), gradgtilde(:), fverif(:), ppformwork(:,:,:) INTEGER, SAVE:: loc_zspan TYPE(mumps_mat), SAVE :: femat !< Finite Element Method matrix for the full domain TYPE(mumps_mat), SAVE :: reduccedmat !< Finite Element Method matrix in the redduced web-spline sub-space TYPE(mumps_mat), SAVE :: fematmpi !< Finite Element Method matrix prepared for mpi parallelism INTEGER :: nbmoments = 10 !< number of moments to be calculated and stored INTEGER(kind=omp_lock_kind), Allocatable:: mu_lock(:) !< Stores the lock for fields parallelism CONTAINS SUBROUTINE mag_init USE basic, ONLY: nr, nz,lu_in,cstep USE magnet ALLOCATE (Br((nr + 1)*(nz + 1)), Bz((nr + 1)*(nz + 1))) ALLOCATE (Athet((nr + 1)*(nz + 1))) ! Calculate magnetic field mirror components in grid points (Davidson analytical formula employed) ! or load it from magnetfile if present call timera(0, "mag_init") CALL magnet_init(lu_in,cstep) call timera(1, "mag_init") end subroutine mag_init !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Set-up the necessary variables for solving Poisson and computes the magnetic field on the grid ! !--------------------------------------------------------------------------- SUBROUTINE fields_init USE basic, ONLY: pot, nlperiod, nrank, rhs, volume, rgrid, Dimensionality,nthet,thetgrid USE bsplines USE geometry USE mumps_bsplines USE mpihelper INTEGER :: nrz(2), i, d2, k1, n1 select case(Dimensionality) case (1) ! Auxiliary vectors ALLOCATE(vec1((nz+1)*(nr+1)),vec2((nr+1)*(nz+1))) DO i=0,nr vec1(i*(nz+1)+1:(i+1)*(nz+1))=zgrid!(0:nz) vec2(i*(nz+1)+1:(i+1)*(nz+1))=rgrid(i) END DO ! Set up 2d spline splrz used in the FEM CALL set_spline(femorder, ngauss, zgrid, rgrid, splrz, nlppform=nlppform, period=nlperiod) ! Set up 2d spline splrz_ext used in the FEM to calculate the external electric field and potential CALL set_spline(femorder, ngauss, zgrid, rgrid, splrz_ext, nlppform=nlppform, period=nlperiod) !Allocate the work buffer to calculate the ppform d2 = splrz%sp2%dim k1 = splrz%sp1%order n1 = splrz%sp1%nints case(2) ALLOCATE(vec1((nthet+1)*(nr+1)),vec2((nr+1)*(nthet+1))) DO i=0,nthet vec1(i*(nr+1)+1:(i+1)*(nr+1))=rgrid!(0:nz) vec2(i*(nr+1)+1:(i+1)*(nr+1))=thetgrid(i) END DO ! Set up 2d spline splrz used in the FEM CALL set_spline(femorder(2:3), ngauss(2:3), rgrid, thetgrid, splrthet, nlppform=nlppform, period=nlperiod(2:3)) ! Set up 2d spline splrz_ext used in the FEM to calculate the external electric field and potential CALL set_spline(femorder(2:3), ngauss(2:3), rgrid, thetgrid, splrthet_ext, nlppform=nlppform, period=nlperiod(2:3)) !Allocate the work buffer to calculate the ppform d2 = splrz%sp2%dim k1 = splrz%sp1%order n1 = splrz%sp1%nints end select ALLOCATE(ppformwork(d2,k1,n1)) ! Calculate dimension of splines nrz(1) = nz nrz(2) = nr CALL get_dim(splrz, nrank, nrz, femorder) ! Allocate necessary variables ALLOCATE (matcoef(nrank(1), nrank(2))) ALLOCATE (pot((nr + 1)*(nz + 1))) ALLOCATE (potxt((nr + 1)*(nz + 1))) ALLOCATE (Erxt((nr + 1)*(nz + 1))) ALLOCATE (Ezxt((nr + 1)*(nz + 1))) ALLOCATE (rhs(nrank(1)*nrank(2))) ALLOCATE (gradgtilde(nrank(1)*nrank(2))) gradgtilde = 0 ALLOCATE (phi_spline(nrank(1)*nrank(2))) ALLOCATE (volume(nrank(1)*nrank(2))) volume = 0 ALLOCATE (Er((nr + 1)*(nz + 1)), Ez((nr + 1)*(nz + 1))) ALLOCATE (mu_lock(nrank(1)*nrank(2))) do i = 1, nrank(1)*nrank(2) call omp_init_lock(mu_lock(i)) end do end SUBROUTINE fields_init !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Set-up the geometry definition and read it from the standard input !> Precomputes the LHS matrix to solve Poisson abd the RHS effect of the dirichlet boundaries ! !--------------------------------------------------------------------------- SUBROUTINE fields_start USE geometry USE basic, ONLY: nrank implicit none INTEGER:: i,j, k, ierr DOUBLE PRECISION:: val ! set up the geometry module for setting up non-conforming boundary conditions call timera(0, "geom_init") call geom_init(splrz, vec1, vec2) call timera(1, "geom_init") ! Initialisation of FEM matrix CALL init(nrank(1)*nrank(2), 2, femat) ! Calculate and factorise FEM matrix (depends only on mesh) CALL fematrixrz(femat) If (walltype .lt. 0) then allocate (fverif(nrank(1)*nrank(2))) fverif = 0 end if ! Compute the volume of the splines and gtilde for solving E using web-splines CALL comp_volume !$OMP PARALLEL Call comp_gradgtilde !$OMP END PARALLEL if (nlweb) then ! Calculate reduced matrix for use of web splines call timera(0, "reduce femat") call Reducematrix(femat, reduccedmat) call timera(1, "reduce femat") - call factor(reduccedmat,nlmetis=.true.) + call factor(reduccedmat) !call init(reduccedmat%rank,2,fematmpi,comm_in=MPI_COMM_WORLD) !do i=1,reduccedmat%rank ! do k=reduccedmat%irow(i),reduccedmat%irow(i+1)-1 ! j=reduccedmat%cols(k) ! call putele(fematmpi, i, j, reduccedmat%val(k)) ! end do !end do else call factor(femat) !call init(femat%rank,2,fematmpi,comm_in=MPI_COMM_WORLD) !do i=1,femat%rank ! do k=femat%irow(i),femat%irow(i+1)-1 ! j=femat%cols(k) ! call putele(fematmpi, i, j, femat%val(k)) ! end do !end do end if !call factor(fematmpi) !WRITE(*,*) "fematmpi is factorised" !WRITE(*,*) "Copy and to_mat worked" !CALL MPI_abort(MPI_COMM_WORLD,-1,ierr) call vacuum_field END SUBROUTINE fields_start !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Recomputes the vacuum electric field ! !--------------------------------------------------------------------------- subroutine vacuum_field Use geometry USE basic, ONLY: pot, rhs implicit none INTEGER:: i, iend ! Computes the externally imposed electric field !$OMP PARALLEL private(i) !$OMP DO SIMD do i=1,nrank(1)*nrank(2) rhs(i)=-gradgtilde(i) !rhs = -gradgtilde if (walltype .lt. 0) rhs(i) = rhs(i) + fverif(i) end do !$OMP END DO SIMD !$OMP END PARALLEL call poisson(splrz_ext) !$OMP PARALLEL private(i,iend) call Update_phi(splrz_ext) !$OMP BARRIER !$OMP DO ! On the root process, compute the electric field for diagnostic purposes DO i=1,size(pot),16 iend=min(size(pot),i+15) potxt(i:iend) = pot(i:iend) erxt(i:iend) = Er(i:iend) Ezxt(i:iend) = Ez(i:iend) END DO !$OMP END DO NOWAIT !$OMP END PARALLEL end subroutine !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Set-up the necessary variables for the communication of moments and rhs grid ! !--------------------------------------------------------------------------- SUBROUTINE fields_comm_init(Zbounds) USE basic, ONLY: nrank USE mpihelper USE omp_lib INTEGER:: Zbounds(0:), nbthreads nbthreads = omp_get_max_threads() Write(*,*) "nbthreads", nbthreads loc_zspan = Zbounds(mpirank + 1) - Zbounds(mpirank) + femorder(1) if (allocated(loc_moments)) deallocate (loc_moments) ALLOCATE (loc_moments(nbmoments, loc_zspan*nrank(2))) if (allocated(loc_rhs)) deallocate (loc_rhs) ALLOCATE (loc_rhs(loc_zspan*nrank(2))) if (allocated(omp_loc_rhs)) deallocate (omp_loc_rhs) ALLOCATE (omp_loc_rhs(nbthreads,loc_zspan*nrank(2))) Write(*,*) "size omp_loc_rhs", size(omp_loc_rhs,1),size(omp_loc_rhs,2) IF (mpisize .gt. 1) THEN CALL init_overlaps(nrank, femorder, Zbounds(mpirank), Zbounds(mpirank + 1), nbmoments) END IF END SUBROUTINE fields_comm_init !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Construct the right hand side vector used in the FEM Poisson solver ! !> @param[in] plist list of the particles type storing the desired specie parameters ! !--------------------------------------------------------------------------- SUBROUTINE rhscon(plist) USE bsplines use mpi USE basic, ONLY: rhs, Zbounds USE beam, ONLY: particles USE mpihelper Use geometry Use omp_lib type(particles), INTENT(INOUT):: plist(:) INTEGER:: i,j,k IF (nlphis) then ! We calculate the self-consistent field !$OMP DO SIMD Do i=1,size(loc_rhs) loc_rhs(i)=0 end do !$OMP END DO SIMD ! Assemble rhs for each specie Do i = 1, size(plist, 1) if (plist(i)%is_field) CALL deposit_charge(plist(i), loc_rhs) END Do !$OMP BARRIER !Communicate the overlaps if(mpisize .gt. 1) call rhs_overlap ! Add gradgtilde !$OMP DO SIMD Do i=0,size(loc_rhs)-1 j=(i)/loc_zspan k=mod(i,loc_zspan) loc_rhs(i+1)=loc_rhs(i+1)-gradgtilde((j)*nrank(1)+(k+Zbounds(mpirank)+1)) end do !$OMP END DO SIMD !add the fverif source for test cases if (walltype .lt. 0)then !$OMP DO Do i=0,size(loc_rhs)-1 j=i/loc_zspan k=mod(i,loc_zspan) loc_rhs(i+1)=loc_rhs(i+1)+fverif((j)*nrank(1)+(k+Zbounds(mpirank)+1)) end do !$OMP END DO end if ! If we are using MPI parallelism, reduce the rhs on the root process IF (mpisize .gt. 1) THEN CALL rhs_gather(rhs) ELSE !in serial copy loc_rhs to rhs !$OMP DO Do i=1,size(loc_rhs) rhs(i)=loc_rhs(i) end do !$OMP END DO END IF ELSE ! We only consider the externally imposed field !$OMP DO Do i=1,size(rhs) rhs(i)=-gradgtilde(i) end do !$OMP END DO END IF END SUBROUTINE rhscon !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculate the 0th 1st and 2nd order moments of the particle p and stores it in moment ! !> @param[in] p the particles type storing the desired specie parameters !> @param[out] moment the 2d array storing the calculated moments ! !--------------------------------------------------------------------------- SUBROUTINE momentsdiag(p) USE bsplines use mpi USE beam, ONLY: particles USE mpihelper Use geometry type(particles), INTENT(INOUT):: p !REAL(kind=db), INTENT(INOUT):: moment(:, :) !$OMP SINGLE loc_moments = 0 ! Reset the moments matrix ! Assemble rhs !$OMP END SINGLE IF (p%Nploc .ne. 0) THEN CALL deposit_moments(p, loc_moments) END IF !$OMP SINGLE if(.not. allocated(p%moments))THEN if(mpirank.eq.0)THEN Allocate(p%moments(nbmoments,nrank(1)*nrank(2))) else Allocate(p%moments(0,0)) end if end if !$OMP END SINGLE ! If we are using MPI parallelism, reduce the rhs on the root process IF (mpisize .gt. 1) THEN CALL moments_gather(p%moments) ELSE !$OMP SINGLE p%moments = loc_moments !$OMP END SINGLE NOWAIT END IF END SUBROUTINE momentsdiag !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Deposit the particles moments (n,v,v^2) from p on the grid ! !> @param[in] p the particles type storing the desired specie parameters !> @param[in] p_loc_moments local tensor used to store the moments of the given specie !--------------------------------------------------------------------------- SUBROUTINE deposit_moments(p, p_loc_moments) USE bsplines use mpi USE basic, ONLY: Zbounds USE beam, ONLY: particles USE mpihelper USE geometry USE omp_lib TYPE(particles), INTENT(IN):: p REAL(kind=db), DIMENSION(:, :), INTENT(INOUT):: p_loc_moments REAL(kind=db), DIMENSION(:, :), Allocatable:: omp_loc_moments INTEGER ::irow, jcol, it, jw, mu, i, k, iend, nbunch INTEGER, DIMENSION(:), ALLOCATABLE::zleft, rleft REAL(kind=db) :: vr, vthet, vz, coeff REAL(kind=db), ALLOCATABLE :: fun(:, :, :), fun2(:, :, :) INTEGER:: num_threads num_threads = omp_get_max_threads() nbunch = p%Nploc/num_threads ! Particle bunch size used when calling basfun nbunch = max(nbunch, 1) ! Particle bunch size used when calling basfun nbunch = min(nbunch, 64) ! Particle bunch size used when calling basfun ! Assemble rhs IF (p%Nploc .gt. 0) THEN !!$OMP PARALLEL DEFAULT(SHARED), PRIVATE(zleft,rleft,jw,it,iend,irow,jcol,mu,k,vr,vz,vthet,coeff,fun,fun2) ALLOCATE (zleft(nbunch), rleft(nbunch)) ALLOCATE (fun(1:femorder(1) + 1, 0:0, nbunch), fun2(1:femorder(2) + 1, 0:0, nbunch)) ! Arrays keeping values of b-splines at gauss node !allocate(omp_loc_moments(size(p_loc_moments,1),size(p_loc_moments,2))) !omp_loc_moments=0 !$OMP DO DO i = 1, p%Nploc, nbunch ! Avoid segmentation fault by accessing non relevant data iend = min(i + nbunch - 1, p%Nploc) k = iend - i + 1 ! Localize the particle !CALL locintv(splrz%sp2, p%R(i:iend), rleft(1:k)) !CALL locintv(splrz%sp1, p%Z(i:iend), zleft(1:k)) rleft(1:k) = p%cellindex(1,i:iend) zleft(1:k) = p%cellindex(3,i:iend) ! Compute the value of the splines at the particles positions CALL basfun(p%pos(3,i:iend), splrz%sp1, fun(:, :, 1:k), zleft(1:k) + 1) CALL basfun(p%pos(1,i:iend), splrz%sp2, fun2(:, :, 1:k), rleft(1:k) + 1) DO k = 1, (iend - i + 1) DO jw = 1, (femorder(2) + 1) DO it = 1, (femorder(1) + 1) irow = zleft(k) + it - Zbounds(mpirank) jcol = rleft(k) + jw mu = irow + (jcol - 1)*(loc_zspan) coeff = p%weight*fun(it, 0, k)*fun2(jw, 0, k) ! Add contribution of particle nbunch to rhs grid point mu vr = 0.5*(p%U(1,i + k - 1)/p%Gamma(i + k - 1) + p%Uold(1,i + k - 1)/p%Gammaold(i + k - 1)) vz = 0.5*(p%U(3,i + k - 1)/p%Gamma(i + k - 1) + p%Uold(3,i + k - 1)/p%Gammaold(i + k - 1)) vthet = 0.5*(p%U(2,i + k - 1)/p%Gamma(i + k - 1) + p%Uold(2,i + k - 1)/p%Gammaold(i + k - 1)) call omp_set_lock(mu_lock(mu)) p_loc_moments(1:10,mu)=p_loc_moments(1:10,mu)+coeff*(/1.0_db,vr,vthet,vz, vr*vr, vr*vthet, vr*vz, vthet**2, vthet*vz, vz**2/) call omp_unset_lock(mu_lock(mu)) END DO END DO END DO END DO !!$OMP END PARALLEL DO !$OMP END DO NOWAIT !Do i=1,size(p_loc_moments,2) ! call omp_set_lock(mu_lock(i)) ! p_loc_moments(:,i)=p_loc_moments(:,i)+omp_loc_moments(:,i) ! call omp_unset_lock(mu_lock(i)) !end do !!$OMP END CRITICAL(loc_moments_reduce) DEALLOCATE (fun, fun2, zleft, rleft) END IF END subroutine deposit_moments !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Deposit the particles charges (q) from p on the grid ! !> @param[in] p the particles type storing the desired specie parameters !> @param[in] p_loc_moments local tensor used to store the moments of the given specie !--------------------------------------------------------------------------- SUBROUTINE deposit_charge(p, p_loc_moments) USE bsplines use mpi USE constants USE basic, ONLY: Zbounds, rnorm, phinorm USE beam, ONLY: particles USE mpihelper USE geometry USE omp_lib TYPE(particles), INTENT(IN):: p REAL(kind=db), DIMENSION(:), INTENT(INOUT):: p_loc_moments INTEGER ::irow, jcol, it, jw, mu, i, k, iend, nbunch INTEGER, DIMENSION(:), ALLOCATABLE::zleft, rleft REAL(kind=db), ALLOCATABLE :: fun(:, :, :), fun2(:, :, :) INTEGER:: num_threads, curr_thread real(kind=db):: contrib, chargecoeff num_threads = omp_get_max_threads() nbunch = p%Nploc/num_threads ! Particle bunch size used when calling basfun nbunch = max(nbunch, 1) ! Particle bunch size used when calling basfun nbunch = min(nbunch, 16) ! Particle bunch size used when calling basfun chargecoeff = p%weight*p%q/(2*pi*eps_0*phinorm*rnorm) ! Normalized charge density simulated by each macro particle ! Assemble rhs IF (p%Nploc .gt. 0) THEN !!!$OMP PARALLEL DEFAULT(SHARED), PRIVATE(i,zleft, rleft, jw, it, iend, irow, jcol, mu, k, fun, fun2, contrib) ALLOCATE (zleft(nbunch), rleft(nbunch)) ALLOCATE (fun(1:femorder(1) + 1, 0:0, nbunch), fun2(1:femorder(2) + 1, 0:0, nbunch)) ! Arrays keeping values of b-splines at gauss node !allocate(omp_loc_moments(size(p_loc_moments))) !omp_loc_moments=0 zleft=0 rleft=0 curr_thread=omp_get_thread_num()+1 !$OMP DO SIMD Do i=1,size(loc_rhs) omp_loc_rhs(:,i)=0 end do !$OMP END DO SIMD !$OMP DO DO i = 1, p%Nploc, nbunch ! Avoid segmentation fault by accessing non relevant data iend = min(i + nbunch - 1, p%Nploc) k = iend - i + 1 ! Localize the particle rleft(1:k) = p%cellindex(1,i:iend) zleft(1:k) = p%cellindex(3,i:iend) ! Compute the value of the splines at the particles positions CALL basfun(p%pos(3,i:iend), splrz%sp1, fun, zleft(1:k) + 1) CALL basfun(p%pos(1,i:iend), splrz%sp2, fun2, rleft(1:k) + 1) !CALL geom_weight(p%Z(i:iend),p%R(i:iend),wgeom) DO k = 1, (iend - i + 1) DO jw = 1, (femorder(2) + 1) DO it = 1, (femorder(1) + 1) irow = zleft(k) + it - Zbounds(mpirank) jcol = rleft(k) + jw mu = irow + (jcol - 1)*(loc_zspan) ! Add contribution of particle k to rhs grid point mu contrib = fun(it, 0, k)*fun2(jw, 0, k)*p%geomweight(0,i + k - 1)*chargecoeff !!$OMP ATOMIC UPDATE omp_loc_rhs(curr_thread,mu) = omp_loc_rhs(curr_thread,mu) + contrib !!$OMP END ATOMIC END DO END DO END DO END DO !$OMP END DO DEALLOCATE (fun, fun2, zleft, rleft) !$OMP DO SIMD Do i=1,size(p_loc_moments) Do k=1,num_threads p_loc_moments(i)=p_loc_moments(i)+omp_loc_rhs(k,i) end do end do !$OMP END DO SIMD END IF END subroutine deposit_charge !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Do the communication of the local moment matrices between mpi workers for the overlap grid points !> ! !--------------------------------------------------------------------------- SUBROUTINE rhs_overlap USE mpihelper USE Basic, ONLY: Zbounds, mpirank, leftproc, rightproc INTEGER:: ierr, i, j !$OMP MASTER !WRITE(*,*) mpirank, "wE communicate overlap rhs" CALL rhsoverlapcomm(mpirank, leftproc, rightproc, loc_rhs, nrank, femorder, loc_zspan - femorder(1)) !$OMP END MASTER !$OMP BARRIER IF (mpirank .gt. 0) THEN !$OMP DO SIMD collapse(2) DO j = 1, femorder(1) DO i = 1, nrank(2) loc_rhs((i - 1)*loc_zspan + j) = loc_rhs((i - 1)*loc_zspan + j)& & + rhsoverlap_buffer(nrank(2)*(j - 1) + i) END DO END DO !$OMP END DO SIMD END IF !$OMP BARRIER END SUBROUTINE rhs_overlap !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Do the communication of the local moment matrices between mpi workers to reduce the result on the host ! !--------------------------------------------------------------------------- SUBROUTINE rhs_gather(rhs) USE mpihelper USE Basic, ONLY: Zbounds, mpirank, leftproc, rightproc REAL(kind=db), DIMENSION(:), INTENT(INOUT):: rhs INTEGER:: ierr, i, j INTEGER:: displs(mpisize), counts(mpisize) INTEGER:: overlap_type INTEGER:: rcvoverlap_type displs = Zbounds(0:mpisize - 1) counts = Zbounds(1:mpisize) - Zbounds(0:mpisize - 1) counts(mpisize) = counts(mpisize) + femorder(1) ! Set communication vector type overlap_type = rhsoverlap_type rcvoverlap_type = rcvrhsoverlap_type !$OMP MASTER IF (mpirank .eq. 0) THEN rhs = 0 END IF CALL MPI_GATHERV(loc_rhs, counts(mpirank + 1), rhsoverlap_type, & & rhs, counts, displs, rcvrhsoverlap_type, 0, MPI_COMM_WORLD, ierr) !$OMP END MASTER !$OMP BARRIER END SUBROUTINE rhs_gather !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Do the communication of the local moment matrices between mpi workers for the overlap grid points !> and reduce the result on the host ! !--------------------------------------------------------------------------- SUBROUTINE moments_gather(moment) USE mpihelper USE Basic, ONLY: Zbounds, mpirank, leftproc, rightproc REAL(kind=db), DIMENSION(:, :), INTENT(INOUT):: moment INTEGER:: ierr, i, j INTEGER:: displs(mpisize), counts(mpisize) displs = Zbounds(0:mpisize - 1) counts = Zbounds(1:mpisize) - Zbounds(0:mpisize - 1) counts(mpisize) = counts(mpisize) + femorder(1) !$OMP MASTER CALL momentsoverlapcomm(mpirank, leftproc, rightproc, loc_moments, nrank, femorder, loc_zspan - femorder(1)) !$OMP END MASTER !$OMP BARRIER IF (mpirank .gt. 0) THEN !!$OMP PARALLEL DO SIMD DEFAULT(SHARED) private(i) !$OMP DO SIMD collapse(2) DO j = 1, femorder(1) DO i = 1, nrank(2) loc_moments(1:nbmoments, (i - 1)*loc_zspan + j) = loc_moments(1:nbmoments, (i - 1)*loc_zspan + j)& & + momentsoverlap_buffer(nbmoments*(nrank(2)*(j - 1) + i - 1) + 1:nbmoments*(nrank(2)*(j - 1) + i)) END DO END DO !$OMP END DO SIMD END IF !$OMP MASTER ! Set communication vector type IF (mpirank .eq. 0) THEN moment = 0 END IF CALL MPI_GATHERV(loc_moments, counts(mpirank + 1), momentsoverlap_type, & & moment, counts, displs, rcvmomentsoverlap_type, 0, MPI_COMM_WORLD, ierr) !$OMP END MASTER !$OMP BARRIER END SUBROUTINE moments_gather !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Solves Poisson equation using FEM. Distributes the result on all MPI workers. ! !--------------------------------------------------------------------------- SUBROUTINE poisson(splinevar) USE basic, ONLY: rhs, nrank, pot, nlend USE bsplines, ONLY: spline2d, gridval USE mumps_bsplines, ONLY: bsolve, vmx USE futils Use geometry type(spline2d):: splinevar INTEGER:: ierr, i, j, iend real(kind=db), allocatable::reducedrhs(:) real(kind=db), allocatable:: reducedsol(:), tempcol(:) allocate (reducedrhs(nrank(1)*nrank(2))) allocate (reducedsol(nbreducedspline)) allocate (tempcol(nrank(1)*nrank(2))) !reduccedmat%mumps_par%ICNTL(11)=1 if (nlweb .and. mpirank.eq.0) then ! we use the web-spline reduction for stability reducedrhs = vmx(etilde, rhs) Call bsolve(reduccedmat, reducedrhs(1:nbreducedspline), reducedsol,10) else if(mpirank.eq.0) then CALL bsolve(femat, rhs, phi_spline) end if !$OMP MASTER if (nlweb) then ! we use the web-spline reduction for stability CALL MPI_Bcast(reducedsol, nbreducedspline, db_type, 0, MPI_COMM_WORLD, ierr) else CALL MPI_Bcast(phi_spline, nrank(1)*nrank(2), db_type, 0, MPI_COMM_WORLD, ierr) end if !$OMP END MASTER !$OMP BARRIER if (nlweb) then ! we use the web-spline reduction for stability tempcol = 0 tempcol(1:nbreducedspline) = reducedsol phi_spline = vmx(etildet, tempcol) end if END SUBROUTINE poisson SUBROUTINE poisson_mpi(splinevar) USE basic, ONLY: rhs, nrank, pot, nlend USE bsplines, ONLY: spline2d, gridval USE mumps_bsplines, ONLY: bsolve, vmx USE futils Use geometry type(spline2d):: splinevar INTEGER:: ierr, i, j, iend real(kind=db), allocatable::reducedrhs(:) real(kind=db), allocatable:: reducedsol(:), tempcol(:) allocate (reducedrhs(nrank(1)*nrank(2))) allocate (reducedsol(nbreducedspline)) allocate (tempcol(nrank(1)*nrank(2))) if (nlweb) then ! we use the web-spline reduction for stability reducedrhs = vmx(etilde, rhs) Call bsolve(fematmpi, reducedrhs(1:nbreducedspline), reducedsol) tempcol = 0 tempcol(1:nbreducedspline) = reducedsol !phi_spline = 0 phi_spline = vmx(etildet, tempcol) else CALL bsolve(fematmpi, rhs, phi_spline) !CALL MPI_Bcast(phi_spline, nrank(1)*nrank(2), db_type, 0, MPI_COMM_WORLD, ierr) end if END SUBROUTINE poisson_mpi !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Updates the splinevar variable with the new phi coefficients and calculates !> Phi Er and Ez on the grid ! !--------------------------------------------------------------------------- SUBROUTINE Update_phi(splinevar) USE basic, ONLY: rhs, nrank, pot, nlend USE bsplines, ONLY: spline2d, gridval USE mumps_bsplines, ONLY: bsolve, vmx USE futils Use geometry type(spline2d):: splinevar INTEGER:: ierr, i, j, iend !$OMP DO SIMD collapse(2) Do j=1,nrank(2) Do i=1,nrank(1) matcoef(i,j) = phi_spline((j-1)*nrank(1)+i) END DO END DO !$OMP END DO SIMD ! update the ppform coefficients CALL updt_ppform2d(splinevar, matcoef) !$OMP BARRIER IF (mpirank .eq. 0 .and. (modulo(step, it2d) .eq. 0 .or. nlend)) THEN !$OMP DO ! On the root process, compute the electric field for diagnostic purposes DO i=1,size(pot),16 iend=min(size(pot),i+15) CALL gridval(splinevar, vec1(i:iend), vec2(i:iend), pot(i:iend), (/0, 0/)) CALL gridval(splinevar, vec1(i:iend), vec2(i:iend), Ez(i:iend), (/1, 0/)) CALL gridval(splinevar, vec1(i:iend), vec2(i:iend), Er(i:iend), (/0, 1/)) Ez(i:iend) = -pot(i:iend)*gridwdir(1,i:iend) - Ez(i:iend)*gridwdir(0,i:iend) - gtilde(1,i:iend) Er(i:iend) = -pot(i:iend)*gridwdir(2,i:iend) - Er(i:iend)*gridwdir(0,i:iend) - gtilde(2,i:iend) pot(i:iend) = pot(i:iend)*gridwdir(0,i:iend) + gtilde(0,i:iend) END DO !$OMP END DO NOWAIT END IF END SUBROUTINE Update_phi !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Computes the electric fields and potential at the particles position for particles !> between positions nstart and nend in the list ! !> @param[in] p the particles type storing the desired specie parameters !> @param[in] nstart starting index for the particle list !> @param[in] nend ending index for the particle list !--------------------------------------------------------------------------- SUBROUTINE EFieldscompatparts(p, nstart, nend) Use beam, ONLY: particles Use geometry Use splinebound TYPE(particles), INTENT(INOUT):: p INTEGER, OPTIONAL::nstart, nend INTEGER:: i, iend, nst, nnd INTEGER:: nbunch INTEGER:: num_threads Real(kind=db), ALLOCATABLE:: erext(:), ezext(:), gtildeloc(:, :) if (.not. present(nstart)) nst = 1 if (.not. present(nend)) nnd = p%Nploc !num_threads = omp_get_max_threads() !nbunch = (nnd - nst + 1)/num_threads ! Particle bunch size used when calling basfun !nbunch = max(nbunch, 1) ! Particle bunch size used when calling basfun nbunch = 64 ! Particle bunch size used when calling basfun Allocate (erext(nbunch), ezext(nbunch), gtildeloc(0:2,0:nbunch - 1)) ! Evaluate the electric potential and field at the particles position !$OMP DO SIMD DO i = nst, nnd, nbunch ! Avoid segmentation fault by accessing non relevant data iend = min(i + nbunch - 1, nnd) CALL speval(splrz, p%pos(3,i:iend), p%pos(1,i:iend),p%cellindex(3,i:iend),p%cellindex(1,i:iend), p%pot(i:iend), p%E(2,i:iend), p%E(1,i:iend)) CALL speval(splrz_ext, p%pos(3,i:iend), p%pos(1,i:iend),p%cellindex(3,i:iend),p%cellindex(1,i:iend), p%potxt(i:iend)) Call total_gtilde(p%pos(3,i:iend), p%pos(1,i:iend), gtildeloc(:,0:iend - i),p%geomweight(:,i:iend)) p%E(2,i:iend) = -p%E(2,i:iend)*p%geomweight(0,i:iend) - p%pot(i:iend)*p%geomweight(1,i:iend) - gtildeloc(1,0:iend - i) p%E(1,i:iend) = -p%E(1,i:iend)*p%geomweight(0,i:iend) - p%pot(i:iend)*p%geomweight(2,i:iend) - gtildeloc(2,0:iend - i) p%pot(i:iend) = p%geomweight(0,i:iend)*p%pot(i:iend) + gtildeloc(0,0:iend - i) p%potxt(i:iend) = p%geomweight(0,i:iend)*p%potxt(i:iend) + gtildeloc(0,0:iend - i) END DO !$OMP END DO SIMD NOWAIT END SUBROUTINE EFieldscompatparts !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Constucts the FEM matrix using bsplines initialized in fields_init !--------------------------------------------------------------------------- SUBROUTINE fematrixrz(mat) USE bsplines USE geometry USE omp_lib USE sparse type(mumps_mat):: mat REAL(kind=db), ALLOCATABLE :: xgauss(:, :), wgauss(:), wgeom(:, :) INTEGER, ALLOCATABLE :: f(:, :), aux(:) REAL(kind=db), ALLOCATABLE :: coefs(:) REAL(kind=db), ALLOCATABLE :: fun(:, :, :), fun2(:, :, :) REAL(kind=db) :: contrib INTEGER, ALLOCATABLE :: idert(:, :), iderw(:, :), iderg(:, :) integer,allocatable:: iid(:),jid(:) INTEGER :: i, j, jt, iw, irow, jcol, mu, igauss, iterm, irow2, jcol2, mu2, kterms, gausssize kterms=8 If (allocated(fun)) deallocate (fun) If (allocated(fun2)) deallocate (fun2) ALLOCATE (fun(1:femorder(1) + 1, 0:1,3*ngauss(1)*ngauss(2)), fun2(1:femorder(2) + 1, 0:1,3*ngauss(1)*ngauss(2))) If (allocated(wgeom)) deallocate (wgeom) ALLOCATE (wgeom(0:2,3*ngauss(1)*ngauss(2)))!Arrays keeping values of b-splines at gauss node !ALLOCATE(xgauss(ngauss(1)*ngauss(2),2), wgauss(ngauss(1)*ngauss(2)),zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2))) !Gaussian nodes and weights arrays ALLOCATE (f((femorder(1) + 1)*(femorder(2) + 1), 2), aux(femorder(1) + 1)) !Auxiliary arrays ordering bsplines ALLOCATE (idert(kterms, 2), iderw(kterms, 2), coefs(kterms), iderg(kterms, 2)) ALLOCATE (iid(3*ngauss(1)*ngauss(2)), jid(3*ngauss(1)*ngauss(2))) !Pointers on the order of derivatives call timera(0, "fematrix") ! Constuction of auxiliary array ordering bsplines in given interval DO i = 1, (femorder(1) + 1) aux(i) = i END DO DO i = 1, (femorder(2) + 1) f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 1) = aux f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 2) = i END DO CALL coefeq(splrz%sp2%knots(0:1), idert, iderw, iderg, coefs, kterms) ! Assemble FEM matrix !$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(j,i,xgauss,wgauss,gausssize,wgeom, igauss,jt,irow,jcol, mu, iw, irow2,jcol2, mu2, contrib, fun, fun2,iid,jid), collapse(2) DO j = 1, nr ! Loop on r position DO i = 1, nz ! Loop on z position !! Computation of gauss weight and position in r and z direction for gaussian integration Call calc_gauss(splrz, ngauss, i, j, xgauss, wgauss, gausssize) iid=i jid=j if (gausssize .gt. 1) then !If (allocated(wgeom)) deallocate (wgeom) !ALLOCATE (wgeom(0:2,gausssize)) CALL geom_weight(xgauss(1:gausssize, 1), xgauss(1:gausssize, 2), wgeom(:,1:gausssize)) CALL basfun(xgauss(1:gausssize, 1), splrz%sp1, fun(:,:,1:gausssize), iid(1:gausssize)) CALL basfun(xgauss(1:gausssize, 2), splrz%sp2, fun2(:,:,1:gausssize), jid(1:gausssize)) End if DO jt = 1, (1 + femorder(1))*(femorder(2) + 1) irow = i + f(jt, 1) - 1; jcol = j + f(jt, 2) - 1 mu = irow + (jcol - 1)*nrank(1) DO iw = 1, (1 + femorder(1))*(femorder(2) + 1) irow2 = i + f(iw, 1) - 1; jcol2 = j + f(iw, 2) - 1 mu2 = irow2 + (jcol2 - 1)*nrank(1) contrib=0.0_db DO igauss = 1, gausssize ! Loop on gaussian weights and positions DO iterm = 1, kterms ! Loop on the two integration dimensions contrib = contrib+wgeom(iderg(iterm, 1),igauss)*wgeom(iderg(iterm, 2),igauss)* & & fun(f(jt, 1), idert(iterm, 1),igauss)*fun(f(iw, 1), idert(iterm, 2),igauss)* & & fun2(f(jt, 2), iderw(iterm, 1),igauss)*fun2(f(iw, 2), iderw(iterm, 2),igauss)* & & wgauss(igauss)*xgauss(igauss, 2) END DO end do call omp_set_lock(mu_lock(mu)) CALL updt_sploc(mat%mat%row(mu), mu2, contrib) call omp_unset_lock(mu_lock(mu)) END DO END DO END DO END DO !$OMP End parallel do DEALLOCATE (f, aux) DEALLOCATE (idert, iderw, coefs, fun, fun2) call timera(1, "fematrix") END SUBROUTINE fematrixrz !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Constucts the FEM matrix using bsplines initialized in fields_init !--------------------------------------------------------------------------- SUBROUTINE fematrixrtheta(mat) USE bsplines USE geometry USE omp_lib USE sparse type(mumps_mat):: mat REAL(kind=db), ALLOCATABLE :: xgauss(:, :), wgauss(:), wgeom(:, :) INTEGER, ALLOCATABLE :: f(:, :), aux(:) REAL(kind=db), ALLOCATABLE :: coefs(:) REAL(kind=db), ALLOCATABLE :: fun(:, :, :), fun2(:, :, :) REAL(kind=db) :: contrib INTEGER, ALLOCATABLE :: idert(:, :), iderw(:, :), iderg(:, :) integer,allocatable:: iid(:),jid(:) INTEGER :: i, j, jt, iw, irow, jcol, mu, igauss, iterm, irow2, jcol2, mu2, kterms, gausssize kterms=8 If (allocated(fun)) deallocate (fun) If (allocated(fun2)) deallocate (fun2) ALLOCATE (fun(1:femorder(1) + 1, 0:1,3*ngauss(1)*ngauss(2)), fun2(1:femorder(2) + 1, 0:1,3*ngauss(1)*ngauss(2))) If (allocated(wgeom)) deallocate (wgeom) ALLOCATE (wgeom(0:2,3*ngauss(1)*ngauss(2)))!Arrays keeping values of b-splines at gauss node !ALLOCATE(xgauss(ngauss(1)*ngauss(2),2), wgauss(ngauss(1)*ngauss(2)),zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2))) !Gaussian nodes and weights arrays ALLOCATE (f((femorder(1) + 1)*(femorder(2) + 1), 2), aux(femorder(1) + 1)) !Auxiliary arrays ordering bsplines ALLOCATE (idert(kterms, 2), iderw(kterms, 2), coefs(kterms), iderg(kterms, 2)) ALLOCATE (iid(3*ngauss(1)*ngauss(2)), jid(3*ngauss(1)*ngauss(2))) !Pointers on the order of derivatives call timera(0, "fematrix") ! Constuction of auxiliary array ordering bsplines in given interval DO i = 1, (femorder(1) + 1) aux(i) = i END DO DO i = 1, (femorder(2) + 1) f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 1) = aux f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 2) = i END DO CALL coefeq(splrz%sp2%knots(0:1), idert, iderw, iderg, coefs, kterms) ! Assemble FEM matrix !$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(j,i,xgauss,wgauss,gausssize,wgeom, igauss,jt,irow,jcol, mu, iw, irow2,jcol2, mu2, contrib, fun, fun2,iid,jid), collapse(2) DO j = 1, nr ! Loop on r position DO i = 1, nz ! Loop on z position !! Computation of gauss weight and position in r and z direction for gaussian integration Call calc_gauss(splrz, ngauss, i, j, xgauss, wgauss, gausssize) iid=i jid=j if (gausssize .gt. 1) then !If (allocated(wgeom)) deallocate (wgeom) !ALLOCATE (wgeom(0:2,gausssize)) CALL geom_weight(xgauss(1:gausssize, 1), xgauss(1:gausssize, 2), wgeom(:,1:gausssize)) CALL basfun(xgauss(1:gausssize, 1), splrz%sp1, fun(:,:,1:gausssize), iid(1:gausssize)) CALL basfun(xgauss(1:gausssize, 2), splrz%sp2, fun2(:,:,1:gausssize), jid(1:gausssize)) End if DO jt = 1, (1 + femorder(1))*(femorder(2) + 1) irow = i + f(jt, 1) - 1; jcol = j + f(jt, 2) - 1 mu = irow + (jcol - 1)*nrank(1) DO iw = 1, (1 + femorder(1))*(femorder(2) + 1) irow2 = i + f(iw, 1) - 1; jcol2 = j + f(iw, 2) - 1 mu2 = irow2 + (jcol2 - 1)*nrank(1) contrib=0.0_db DO igauss = 1, gausssize ! Loop on gaussian weights and positions DO iterm = 1, kterms ! Loop on the two integration dimensions contrib = contrib+wgeom(iderg(iterm, 1),igauss)*wgeom(iderg(iterm, 2),igauss)* & & fun(f(jt, 1), idert(iterm, 1),igauss)*fun(f(iw, 1), idert(iterm, 2),igauss)* & & fun2(f(jt, 2), iderw(iterm, 1),igauss)*fun2(f(iw, 2), iderw(iterm, 2),igauss)* & & wgauss(igauss)*xgauss(igauss, 1) END DO end do call omp_set_lock(mu_lock(mu)) CALL updt_sploc(mat%mat%row(mu), mu2, contrib) call omp_unset_lock(mu_lock(mu)) END DO END DO END DO END DO !$OMP End parallel do DEALLOCATE (f, aux) DEALLOCATE (idert, iderw, coefs, fun, fun2) call timera(1, "fematrix") END SUBROUTINE fematrixrtheta !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Computes the volume of the splines cells needed to display the density in post-processing !--------------------------------------------------------------------------- SUBROUTINE comp_volume USE bsplines USE geometry USE basic, ONLY: Volume REAL(kind=db), ALLOCATABLE :: xgauss(:, :), wgauss(:), wgeom(:, :) INTEGER, ALLOCATABLE :: f(:, :), aux(:) REAL(kind=db), ALLOCATABLE :: coefs(:) REAL(kind=db), ALLOCATABLE :: fun(:, :), fun2(:, :), ftestpt(:, :) Integer, ALLOCATABLE, Dimension(:) :: idg, idt, idp, idw INTEGER :: i, j, jt, irow, jcol, mu, igauss, gausssize, iterm, nterms Real(kind=db)::newcontrib call timera(0, "comp_volume") ALLOCATE (fun(1:femorder(1) + 1, 0:1), fun2(1:femorder(2) + 1, 0:1))!Arrays keeping values of b-splines at gauss node !ALLOCATE(xgauss(ngauss(1)*ngauss(2),2), wgauss(ngauss(1)*ngauss(2)),zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2))) !Gaussian nodes and weights arrays ALLOCATE (f((femorder(1) + 1)*(femorder(2) + 1), 2), aux(femorder(1) + 1)) !Auxiliary arrays ordering bsplines nterms = 4 Allocate (idg(nterms), idt(nterms), idw(nterms), idp(nterms), coefs(nterms)) ! Constuction of auxiliary array ordering bsplines in given interval DO i = 1, (femorder(1) + 1) aux(i) = i END DO DO i = 1, (femorder(2) + 1) f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 1) = aux f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 2) = i END DO volume = 0 if (walltype .lt. 0) fverif = 0 ! Assemble Volume matrix !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(j,i,xgauss,wgauss,gausssize,wgeom, igauss, ftestpt, iterm,jt,irow,jcol, mu, idw, idt, idg, idp, coefs, fun, fun2, newcontrib), collapse(2) DO j = 1, nr ! Loop on r position DO i = 1, nz ! Loop on z position ! Computation of gauss weight and position in r and z direction for gaussian integration Call calc_gauss(splrz, ngauss, i, j, xgauss, wgauss, gausssize) If (allocated(wgeom)) deallocate (wgeom) if (gausssize .gt. 0) then ALLOCATE (wgeom(0:2,size(xgauss, 1))) CALL geom_weight(xgauss(:, 1), xgauss(:, 2), wgeom) End if if (walltype .lt. 0) then If (allocated(ftestpt)) deallocate (ftestpt) ALLOCATE (ftestpt(0:0,size(xgauss, 1))) CALL ftest(xgauss(:, 1), xgauss(:, 2), ftestpt) end if DO igauss = 1, gausssize ! Loop on gaussian weights and positions CALL basfun(xgauss(igauss, 1), splrz%sp1, fun, i) CALL basfun(xgauss(igauss, 2), splrz%sp2, fun2, j) !CALL coefeqext(xgauss(igauss, :), idt, idw, idg, idp, coefs) DO jt = 1, (1 + femorder(1))*(femorder(2) + 1) irow = i + f(jt, 1) - 1; jcol = j + f(jt, 2) - 1 mu = irow + (jcol - 1)*nrank(1) newcontrib = 2*pi*fun(f(jt, 1), 0)*fun2(f(jt, 2), 0)*wgauss(igauss)*xgauss(igauss, 2)!*wgeom(igauss,0) !$OMP ATOMIC UPDATE volume(mu) = volume(mu) + newcontrib !$OMP END ATOMIC if (walltype .lt. 0) THEN newcontrib = ftestpt(0,igauss)*fun(f(jt, 1), 0)*fun2(f(jt, 2), 0)& &*wgeom(0,igauss)*wgauss(igauss)*xgauss(igauss, 2) !$OMP ATOMIC UPDATE fverif(mu) = fverif(mu) + newcontrib !$OMP END ATOMIC end if END DO END DO END DO END DO !$OMP END PARALLEL DO !DEALLOCATE(xgauss, wgauss,zg,rg, wzg, wrg) DEALLOCATE (f, aux) DEALLOCATE (fun, fun2) call timera(1, "comp_volume") END SUBROUTINE comp_volume !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Computes the gradient of the gtilde function for the web-spline method needed to correctly apply the dirichlet boundary conditions !--------------------------------------------------------------------------- SUBROUTINE comp_gradgtilde USE bsplines USE geometry REAL(kind=db), ALLOCATABLE :: xgauss(:, :), wgauss(:), wgeom(:, :) INTEGER, ALLOCATABLE :: f(:, :), aux(:) REAL(kind=db), ALLOCATABLE :: coefs(:) REAL(kind=db), ALLOCATABLE :: fun(:, :,:), fun2(:, :,:), gtildeintegr(:, :) Integer, ALLOCATABLE, Dimension(:) :: idg, idt, idp, idw integer,allocatable:: iid(:),jid(:) INTEGER :: i, j, jt, irow, jcol, mu, igauss, gausssize, iterm, nterms Real(kind=db)::newcontrib !call timera(0, "comp_gradgtilde") ALLOCATE (fun(1:femorder(1) + 1, 0:1,3*ngauss(1)*ngauss(2)), fun2(1:femorder(2) + 1, 0:1,3*ngauss(1)*ngauss(2)))!Arrays keeping values of b-splines at gauss node ALLOCATE (wgeom(0:2,3*ngauss(1)*ngauss(2)))!Arrays keeping values of b-splines at gauss node ALLOCATE (f((femorder(1) + 1)*(femorder(2) + 1), 2), aux(femorder(1) + 1)) !Auxiliary arrays ordering bsplines nterms = 4 Allocate (idg(nterms), idt(nterms), idw(nterms), idp(nterms), coefs(nterms)) ALLOCATE (iid(3*ngauss(1)*ngauss(2)), jid(3*ngauss(1)*ngauss(2))) If (allocated(gtildeintegr)) deallocate (gtildeintegr) ALLOCATE (gtildeintegr(0:2,3*ngauss(1)*ngauss(2))) ! Constuction of auxiliary array ordering bsplines in given interval DO i = 1, (femorder(1) + 1) aux(i) = i END DO DO i = 1, (femorder(2) + 1) f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 1) = aux f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 2) = i END DO CALL coefeqext(splrz%sp2%knots(0:1), idt, idw, idg, idp, coefs) !$OMP DO SIMD do j=1,size(gradgtilde) gradgtilde(j) = 0 END DO !$OMP END DO SIMD !$OMP BARRIER ! Assemble gradgtilde matrix !$OMP DO collapse(2) DO j = 1, nr ! Loop on r position DO i = 1, nz ! Loop on z position ! Computation of gauss weight and position in r and z direction for gaussian integration Call calc_gauss(splrz, ngauss, i, j, xgauss, wgauss, gausssize) iid=i jid=j if (gausssize .gt. 1) then !If (allocated(wgeom)) deallocate (wgeom) !ALLOCATE (wgeom(0:2,gausssize)) CALL geom_weight(xgauss(1:gausssize, 1), xgauss(1:gausssize, 2), wgeom(:,1:gausssize)) CALL basfun(xgauss(1:gausssize, 1), splrz%sp1, fun(:,:,1:gausssize), iid(1:gausssize)) CALL basfun(xgauss(1:gausssize, 2), splrz%sp2, fun2(:,:,1:gausssize), jid(1:gausssize)) Call total_gtilde(xgauss(1:gausssize, 1), xgauss(1:gausssize, 2), gtildeintegr(:,1:gausssize),wgeom(:,1:gausssize)) End if DO jt = 1, (1 + femorder(1))*(femorder(2) + 1) irow = i + f(jt, 1) - 1; jcol = j + f(jt, 2) - 1 mu = irow + (jcol - 1)*nrank(1) newcontrib = 0.0_db DO igauss = 1, gausssize ! Loop on gaussian weights and positions Do iterm = 1, nterms newcontrib = newcontrib + wgeom( idg(iterm),igauss)*gtildeintegr( idp(iterm),igauss)* & & fun(f(jt, 1), idt(iterm),igauss)*fun2(f(jt, 2), idw(iterm),igauss)* & & wgauss(igauss)*xgauss(igauss, 2) End do end do !$OMP ATOMIC UPDATE gradgtilde(mu) = gradgtilde(mu) + newcontrib !$OMP END ATOMIC END DO END DO END DO !$OMP END DO !DEALLOCATE(xgauss, wgauss,zg,rg, wzg, wrg) DEALLOCATE (f, aux) DEALLOCATE (fun, fun2) !call timera(1, "comp_gradgtilde") END SUBROUTINE comp_gradgtilde !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Imposes the dirichlet boundary conditions on the FEM matrix for the case where we use regular splines ( not web-splines). !--------------------------------------------------------------------------- SUBROUTINE fe_dirichlet REAL(kind=db), ALLOCATABLE :: arr(:) INTEGER :: i ALLOCATE (arr(nrank(1)*nrank(2))) DO i = 1, nrank(1) IF (rgrid(0) .ne. 0.0_db) THEN arr = 0; arr(i) = 1; CALL putrow(femat, i, arr) END IF arr = 0; arr(nrank(1)*nrank(2) + 1 - i) = 1; CALL putrow(femat, nrank(1)*nrank(2) + 1 - i, arr) END DO DEALLOCATE (arr) END SUBROUTINE fe_dirichlet !________________________________________________________________________________ SUBROUTINE coefeq(x, idt, idw, idg, c, kterms) REAL(kind=db), INTENT(in) :: x(:) INTEGER, INTENT(out) :: idt(:, :), idw(:, :), idg(:, :),kterms REAL(kind=db), INTENT(out) :: c(:) kterms=8 c = x(2) idt(1, 1) = 0 idt(1, 2) = 0 idw(1, 1) = 0 idw(1, 2) = 0 idg(1, 1) = 1 idg(1, 2) = 1 idt(2, 1) = 0 idt(2, 2) = 1 idw(2, 1) = 0 idw(2, 2) = 0 idg(2, 1) = 1 idg(2, 2) = 0 idt(3, 1) = 1 idt(3, 2) = 0 idw(3, 1) = 0 idw(3, 2) = 0 idg(3, 1) = 0 idg(3, 2) = 1 idt(4, 1) = 1 idt(4, 2) = 1 idw(4, 1) = 0 idw(4, 2) = 0 idg(4, 1) = 0 idg(4, 2) = 0 idt(5, 1) = 0 idt(5, 2) = 0 idw(5, 1) = 0 idw(5, 2) = 0 idg(5, 1) = 2 idg(5, 2) = 2 idt(6, 1) = 0 idt(6, 2) = 0 idw(6, 1) = 0 idw(6, 2) = 1 idg(6, 1) = 2 idg(6, 2) = 0 idt(7, 1) = 0 idt(7, 2) = 0 idw(7, 1) = 1 idw(7, 2) = 0 idg(7, 1) = 0 idg(7, 2) = 2 idt(8, 1) = 0 idt(8, 2) = 0 idw(8, 1) = 1 idw(8, 2) = 1 idg(8, 1) = 0 idg(8, 2) = 0 END SUBROUTINE coefeq SUBROUTINE coefeqext(x, idt, idw, idg, idp, c) REAL(kind=db), INTENT(in) :: x(:) INTEGER, INTENT(out) :: idp(:), idt(:), idw(:), idg(:) REAL(kind=db), INTENT(out) :: c(:) c(1) = x(2) idp(1) = 1 idg(1) = 1 idt(1) = 0 idw(1) = 0 c(2) = x(2) idp(2) = 1 idg(2) = 0 idt(2) = 1 idw(2) = 0 c(3) = x(2) idp(3) = 2 idg(3) = 2 idt(3) = 0 idw(3) = 0 c(4) = x(2) idp(4) = 2 idg(4) = 0 idt(4) = 0 idw(4) = 1 END SUBROUTINE coefeqext - !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Free the memory used by the fields module !--------------------------------------------------------------------------- SUBROUTINE clean_fields Use bsplines USE basic, ONLY: rhs INTEGER:: i do i = 1, nrank(1)*nrank(2) call omp_destroy_lock(mu_lock(i)) end do DEALLOCATE (mu_lock) DEALLOCATE (matcoef) DEALLOCATE (pot) DEALLOCATE (rhs) DEALLOCATE (loc_rhs) DEALLOCATE (loc_moments) DEALLOCATE (phi_spline) DEALLOCATE (Br, Bz) DEALLOCATE (Er, Ez) DEALLOCATE (vec1, vec2) Call DESTROY_SP(splrz) Call DESTROY_SP(splrz_ext) END SUBROUTINE clean_fields SUBROUTINE updt_sploc(arow, j, val) ! ! Update element j of row arow or insert it in an increasing "index" ! USE sparse TYPE(sprow), TARGET :: arow INTEGER, INTENT(in) :: j DOUBLE PRECISION, INTENT(in) :: val ! TYPE(elt), TARGET :: pre_root TYPE(elt), POINTER :: t, p ! if(val.eq.0) return pre_root%next => arow%row0 ! pre_root is linked to the head of the list. t => pre_root DO WHILE (ASSOCIATED(t%next)) p => t%next IF (p%index .EQ. j) THEN p%val = p%val + val RETURN END IF IF (p%index .GT. j) EXIT t => t%next END DO ALLOCATE (p) p = elt(j, val, t%next) t%next => p ! arow%nnz = arow%nnz + 1 arow%row0 => pre_root%next ! In case the head is altered END SUBROUTINE updt_sploc SUBROUTINE updt_ppform2d(sp,c) use bsplines TYPE(spline2d), INTENT(inout) :: sp DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: c !DOUBLE PRECISION, ALLOCATABLE :: work(:,:,:) INTEGER:: m,mm INTEGER :: d1, d2, k1, k2, n1, n2 d1 = sp%sp1%dim d2 = sp%sp2%dim k1 = sp%sp1%order k2 = sp%sp2%order n1 = sp%sp1%nints n2 = sp%sp2%nints !$OMP SINGLE IF( ASSOCIATED(sp%bcoefs) ) DEALLOCATE(sp%bcoefs) ALLOCATE(sp%bcoefs(SIZE(c,1),SIZE(c,2))) !$OMP END SINGLE !ALLOCATE(work(d2,k1,n1)) !$OMP DO DO m=1,SIZE(c,2) CALL topp0(sp%sp1, c(:,m), ppformwork(m,:,:)) sp%bcoefs(:,m)=c(:,m) END DO !$OMP END DO NOWAIT !$OMP SINGLE IF( ASSOCIATED(sp%ppform) ) DEALLOCATE(sp%ppform) ALLOCATE(sp%ppform(k1,n1,k2,n2)) !$OMP END SINGLE !$OMP DO DO mm=1,SIZE(ppformwork,3) DO m=1,SIZE(ppformwork,2) CALL topp0(sp%sp2, ppformwork(:,m,mm), sp%ppform(m,mm,:,:)) END DO END DO !$OMP END DO !DEALLOCATE(work) end subroutine updt_ppform2d !=========================================================================== SUBROUTINE topp0(sp, c, ppform) ! ! Compute PPFORM of a fuction defined by the spline SP ! and spline coefficients C(1:d) ! use bsplines TYPE(spline1d), INTENT(in) :: sp DOUBLE PRECISION, INTENT(in) :: c(:) DOUBLE PRECISION, INTENT(out) :: ppform(0:,:) INTEGER :: p, nints, i, j, k ! p = sp%order - 1 nints = sp%nints ! ppform = 0.0d0 DO i=1,nints ! on each knot interval DO j=1,p+1 ! all spline in interval i DO k=0,p ! k_th derivatives ppform(k,i) = ppform(k,i) + sp%val0(k,j,i)*c(j+i-1) END DO END DO END DO ! END SUBROUTINE topp0 !+ END MODULE fields diff --git a/src/ion_induced_mod.f90 b/src/ion_induced_mod.f90 index 8444e67..4e209c9 100644 --- a/src/ion_induced_mod.f90 +++ b/src/ion_induced_mod.f90 @@ -1,452 +1,452 @@ !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: IIEE ! !> @author !> S. Guinchard - EPFL/SPC ! !> Last modif. !> 11/28 2022 ! ! DESCRIPTION: !> Module handling ion induced electron emissions (IIEE) !> following Schou's model (for Kinetic emissions) and Auger neutralisation !> for potential emissions !------------------------------------------------------------------------------ MODULE iiee USE particletypes USE constants USE basic USE materials !#include "mkl_vsl.f90" ! for random # generators using MKL intel library IMPLICIT NONE CONTAINS !--------------------------------------------------------------------------- ! SUBROUTINE ion_induced(pion, losthole, pelec, nblostparts) ! ! ! DESCRIPTION !> function to determine the number of electrons !> to add to a given species as a function of th number !> of lost ions ! !-------------------------------------------------------------------------- SUBROUTINE ion_induced(pion, losthole, pelec, nblostparts) USE geometry TYPE(particles), INTENT(INOUT):: pion, pelec !< ion and electrons parts REAL(KIND = db), DIMENSION(3) :: last_pos !< last position for lost ion (revert push) REAL(KIND = db), DIMENSION(3) :: normal_dir !< normal direction vector (normalised) INTEGER, DIMENSION(pion%Nploc):: losthole !< indices of lost ions INTEGER ::i,j, nblostparts, Nploc, Nploc_old, Nploc_init!< loop indices and #lost particles INTEGER :: parts_size_increase, nbadded INTEGER :: neuttype_id !< neutral gas type_id INTEGER :: material_id !< electrode material type_id INTEGER :: gen_el, kmax !< # of electrons generated, max# possibly gen. elec. REAL(KIND = db) :: lambda !< Poisson param. to gen elec. (yield) REAL(KIND = db) :: kappa, theta, Emax !< gamma distribution parameters REAL(KIND = db) :: Ekin, Eem !< kinetic energy of lost particles (yield param) and of emitted electrons kmax = 12 !> Max num. elec. to be generated (Poisson) kappa = 4.0 !> kappa param. (Gamma) theta = 0.5 !> theta param. (Gamma) Emax = 25 !> Max value for el. (Gamma) neuttype_id = pion%neuttype_id !> temporarily stored in particle type material_id = pion%material_id !> temporarily stored in particle type IF(pelec%Nploc + 2*nblostparts .gt. size(pelec%pos,2)) THEN parts_size_increase=Max(floor(0.1*size(pelec%pos,2)),2*nblostparts) CALL change_parts_allocation(pelec, parts_size_increase) END IF Nploc_init = pelec%Nploc DO i=1,nblostparts Ekin = compute_Ekin( pion%U(:,losthole(i)), pion) - IF (pion%neuttype_id == 1) THEN ! Yield for H2 + IF (pion%neuttype_id .eq. 1) THEN ! Yield for H2 lambda = 2.0*compute_yield(Ekin/2.0, neuttype_id, material_id) ! computed from yield by protons ELSE ! Yield for other neutral lambda = compute_yield(Ekin, neuttype_id, material_id) ! needs to be changed for other gases END IF nbadded = gen_elec(lambda, kmax) Nploc_old = pelec%Nploc pelec%Nploc = pelec%Nploc + nbadded Nploc = pelec%Nploc last_pos = revert_push(pion, losthole(i)) pelec%nbadded = pelec%nbadded+nbadded normal_dir = find_normal(last_pos) DO j=1,nbadded pelec%pos(1,Nploc_old+j) = last_pos(1) pelec%pos(2,Nploc_old+j) = last_pos(2) pelec%pos(3,Nploc_old+j) = last_pos(3) pelec%newindex = pelec%newindex + 1 pelec%partindex(Nploc_old+j) = pelec%newindex - IF(pelec%zero_vel == .false.) THEN + IF(pelec%zero_vel .eqv. .false.) THEN Eem = gen_E_gamma(kappa, theta, Emax) !> generate an energy value following gamma distribution pelec%U(1,Nploc_old+j) = compute_Vnorm(Eem, pelec)* normal_dir(1) !> Vr pelec%U(3,Nploc_old+j) = compute_Vnorm(Eem, pelec)* normal_dir(3) !> Vz pelec%U(2,Nploc_old+j) = compute_Vnorm(Eem, pelec)* normal_dir(2) !> Vthet ELSE pelec%U(1,Nploc_old+j) = 0.0 pelec%U(3,Nploc_old+j) = 0.0 pelec%U(2,Nploc_old+j) = 0.0 END IF END DO END DO if (pelec%Nploc-Nploc_init .ge. 1) then call geom_weight(pelec%pos(3,Nploc_init+1:pelec%Nploc), pelec%pos(1,Nploc_init+1:pelec%Nploc), pelec%geomweight(:,Nploc_init+1:pelec%Nploc)) end if END SUBROUTINE ion_induced !--------------------------------------------------------------------------- ! FUNCTION compute_Ekin(velocity, p) ! ! ! DESCRIPTION !> Computes the kinetic energy in MeV of a particle given its 3-vel. components !-------------------------------------------------------------------------- FUNCTION compute_Ekin(velocity, p) RESULT(Ekin) TYPE(particles), INTENT(INOUT):: p REAL(KIND = db), DIMENSION(3) :: velocity REAL(KIND = db) :: Ekin Ekin = 5E-7 * p%m * vlight**2 /elchar * (velocity(1)**2 + velocity(2)**2 + velocity(3)**2) END FUNCTION compute_Ekin !--------------------------------------------------------------------------- ! FUNCTION compute_Vnorm(Ekin,p) ! ! ! DESCRIPTION !> Computes the normal velocity of an incident electron emitted !> with energy Ekin !-------------------------------------------------------------------------- FUNCTION compute_Vnorm(Ekin, p) RESULT(Vnorm) REAL(KIND = db) :: Ekin, Vnorm !> Ekin of emitted electron, Normal. corres. veloc. TYPE(particles) :: p !> electrons Vnorm = sqrt(2/p%m * Ekin * elchar) / vlight !> * elchar to get the enery in J and Vnorm in m/s END FUNCTION compute_Vnorm !--------------------------------------------------------------------------- ! FUNCTION fin_normal(last_position) ! ! ! DESCRIPTION !> Computes the normal velocity of an incident electron emitted !> with energy Ekin !-------------------------------------------------------------------------- FUNCTION find_normal(last_position) RESULT(normal_dir) USE geometry REAL(KIND = db), DIMENSION(3) :: last_position !> Last pos. to eval. geom. weight at REAL(KIND = db), DIMENSION(3) :: normal_dir !> Normal direction vector (Result) REAL(KIND = db), DIMENSION(3) :: weight !> Geom. weight at last pos. REAL(KIND = db) :: norm !> To normalise normal vect. call geom_weight(last_position(3), last_position(1), weight) norm = sqrt(weight(2)**2 + weight(3)**2) normal_dir(1) = 1/norm * weight(3) !> Normal along r normal_dir(2) = 0.0 !> Normal along theta normal_dir(3) = 1/norm * weight(2) !> Normal along z END FUNCTION find_normal !--------------------------------------------------------------------------- ! FUNCTION revert_push(pion, partid) ! ! ! DESCRIPTION !> reverts Buneman algorithm over one time step !> to obtain ion position right before recapture !> and hence new electron positions ! !-------------------------------------------------------------------------- FUNCTION revert_push(pion, partid) USE basic, ONLY: dt, tnorm REAL(KIND=db), DIMENSION(3):: revert_push TYPE(particles), INTENT(INOUT):: pion !> species: ions INTEGER :: partid !> id of particle to reverse position revert_push(1) = pion%pos(1,partid) - pion%U(1,partid)*dt revert_push(2) = pion%pos(2,partid) -1/pion%pos(1,partid)* pion%U(2,partid)*dt revert_push(3) = pion%pos(3,partid) -pion%U(3,partid)*dt END FUNCTION revert_push !--------------------------------------------------------------------------- ! FUNCTION eval_polynomial(coefficients, valeur) ! ! ! DESCRIPTION !> Evaluate a polynomial at a given point !> with its coefficients provided in an array !> s.t lowest order coeff = 1st element ! !-------------------------------------------------------------------------- REAL(KIND = db) FUNCTION eval_polynomial(coefficients, valeur) REAL(KIND = db), DIMENSION(:) :: coefficients !< polynomial (e.g fitted yield) coeffs REAL(KIND = db) :: valeur !< point where to evaluate polyn INTEGER :: ii eval_polynomial = 0 DO ii=1, size(coefficients) eval_polynomial = eval_polynomial+coefficients(ii)*valeur**(ii-1) END DO END FUNCTION eval_polynomial !--------------------------------------------------------------------------- ! FUNCTION compute_yield(energy, neuttype_id, material_id) ! ! ! DESCRIPTION !> Gives the theoretical value for the electron yield !> as a function of the energy of the incident ion and !> the type of neutral gas ! !-------------------------------------------------------------------------- REAL(KIND = db) FUNCTION compute_yield(energy, neuttype_id, material_id) !add material id asap REAL(KIND = db) :: energy INTEGER :: neuttype_id, material_id REAL(KIND = db) :: Lambda_exp Lambda_exp = 1E-3 SELECT CASE(material_id) CASE(1) !304 stainless steel IF(energy.le. 1E-3 ) THEN SELECT CASE(neuttype_id) CASE(1) compute_yield = eval_polynomial(coefficients_1H_SS, energy) CASE(2) compute_yield = eval_polynomial(coefficients_1He_SS, energy) CASE(3) compute_yield = eval_polynomial(coefficients_1Ne_SS, energy) CASE DEFAULT compute_yield = eval_polynomial(coefficients_1H_SS, energy) END SELECT ELSE IF(energy.gt. 1E-3 .and. energy.le. 1E-2) THEN compute_yield = Lambda_exp * eval_polynomial(coefficients_2_SS,energy) ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN compute_yield = Lambda_exp * eval_polynomial(coefficients_3_SS,energy) ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN compute_yield = Lambda_exp * eval_polynomial(coefficients_4_SS,energy) ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN compute_yield = Lambda_exp * eval_polynomial(coefficients_5_SS,energy) END IF CASE(2) ! Copper IF(energy.le. 1E-3 ) THEN SELECT CASE(neuttype_id) CASE(1) compute_yield = eval_polynomial(coefficients_1H_Cu, energy) CASE(2) compute_yield = eval_polynomial(coefficients_1He_Cu, energy) CASE(3) compute_yield = eval_polynomial(coefficients_1Ne_Cu, energy) CASE DEFAULT compute_yield = eval_polynomial(coefficients_1H_Cu, energy) END SELECT ELSE IF(energy.gt. 1E-3 .and. energy.le. 1E-2) THEN compute_yield = Lambda_exp * eval_polynomial(coefficients_2_Cu,energy) ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN compute_yield = Lambda_exp * eval_polynomial(coefficients_3_Cu,energy) ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN compute_yield = Lambda_exp * eval_polynomial(coefficients_4_Cu,energy) ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN compute_yield = Lambda_exp * eval_polynomial(coefficients_5_Cu,energy) END IF CASE(3) ! Alumium IF(energy.le. 1E-3 ) THEN SELECT CASE(neuttype_id) CASE(1) compute_yield = eval_polynomial(coefficients_1H_Al, energy) CASE(2) compute_yield = eval_polynomial(coefficients_1He_Al, energy) CASE(3) compute_yield = eval_polynomial(coefficients_1Ne_Al, energy) CASE DEFAULT compute_yield = eval_polynomial(coefficients_1H_Al, energy) END SELECT ELSE IF(energy.gt. 1E-3 .and. energy.le. 1E-2) THEN compute_yield = Lambda_exp * eval_polynomial(coefficients_2_Al,energy) ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN compute_yield = Lambda_exp * eval_polynomial(coefficients_3_Al,energy) ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN compute_yield = Lambda_exp * eval_polynomial(coefficients_4_Al,energy) ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN compute_yield = Lambda_exp * eval_polynomial(coefficients_5_Al,energy) END IF END SELECT END FUNCTION compute_yield !--------------------------------------------------------------------------- ! FUNCTION gen_elec(lambda, kmax) ! ! ! DESCRIPTION !> Gives random values distributed !> following a Poisson distrib. of parameter !> lambda = yield(E) for incomin ion energy E !-------------------------------------------------------------------------- INTEGER FUNCTION gen_elec(lambda, kmax) - USE random + USE random_mod REAL(KIND = db) :: lambda !< Lambda parameter for Poisson distribution REAL(KIND = db) :: nb_alea(1:1) !< random number unif. generated in [0,1] INTEGER :: kmax !< max number possible from Poisson REAL(KIND = db) :: CumulPoisson !< Flag to ensure CDF ~ 1 INTEGER :: i, ii !< loop indices REAL(KIND = db), DIMENSION(kmax) :: vect, SumPart !< terms, partial sums for CDF !> Compute probabilities for each int. value and CDF values DO i = 1,kmax vect(i) = exp(-lambda)*lambda**(i-1)/factorial_fun(i-1); SumPart(i) = sum(vect(1:i)); END DO !> CDF expected to sum to 1 CumulPoisson = sum(vect) !> Generate poisson distrib. int. (see Matlab. code for convg.) call random_array(nb_alea,1,ran_index(1),ran_array(:,1)) DO ii = 1,size(SumPart)-1 IF (nb_alea(1) .lt. SumPart(1)) THEN gen_elec = 0 ELSE IF ((SumPart(ii).le.nb_alea(1)) .and. (nb_alea(1) .lt. SumPart(ii+1))) THEN gen_elec = ii END IF END DO ! Below: see Intel oneAPI Math Kernel Library - Fortran ! to optimise running speed when generating random numbers ! TO DO : change if enough time !----------------------------------------------------------- ! ! INTEGER, INTENT(IN) :: method ! TYPE (VSL_STREAM_STATE), INTENT(IN) :: stream ! INTEGER, INTENT(IN) :: n = 1 ! INTEGER, INTENT(IN) :: r ! status = virngpoisson( method, stream, n, r, lambda ) ! gen_elec = r ! !------------------------------------------------------------ END FUNCTION gen_elec !--------------------------------------------------------------------------- ! FUNCTION gen_E_gamma(kappa, theta, Emax) ! ! ! DESCRIPTION !> Gives random values distributed !> following a Gamma distrib. of parameters (kappa, theta) !> in [0, Emax] eV and peaked at E=2eV !-------------------------------------------------------------------------- FUNCTION gen_E_gamma(kappa, theta, Emax) RESULT(E_el) - USE random + USE random_mod USE incomplete_gamma REAL(KIND = db) :: E_el, Emin, Emax REAL(KIND = db) :: kappa, theta !> parameters to shape Gamma_distr REAL(KIND = db) :: nb_alea(1:1) REAL(KIND = db), DIMENSION(1000) :: Einc, CDF, Diff INTEGER :: ii !> loop index INTEGER :: ifault INTEGER :: posid Emin = 0.00 Emax = 25.0 kappa = 4.0 theta = 0.5 DO ii= 1,size(Einc) Einc(ii) = ( Emin + (ii-1)*(Emax-Emin)/(size(Einc)-1) ) CDF(ii) = gamain(Einc(ii)/theta, kappa, ifault) END DO !> Generate Gamma distrib. int. (see Matlab. code for convg.) call random_array(nb_alea,1,ran_index(1),ran_array(:,1)) Diff = abs(nb_alea(1) - CDF) posid = minloc(Diff, dim = 1) E_el = Einc(posid) END FUNCTION gen_E_gamma !--------------------------------------------------------------------------- ! FUNCTION factorial_fun(n) ! ! ! DESCRIPTION !> Gives the factorial of an integer !-------------------------------------------------------------------------- INTEGER FUNCTION factorial_fun(n) INTEGER :: ii INTEGER :: n factorial_fun = 1 IF (n .lt. 0) THEN RETURN - ELSE IF (n == 0) THEN + ELSE IF (n .eq. 0) THEN factorial_fun = 1 ELSE DO ii=1,n factorial_fun = factorial_fun*ii END DO END IF END FUNCTION factorial_fun END MODULE iiee diff --git a/src/magnet_mod.f90 b/src/magnet_mod.f90 index ab92b54..db44d74 100644 --- a/src/magnet_mod.f90 +++ b/src/magnet_mod.f90 @@ -1,399 +1,399 @@ module magnet use constants IMPLICIT NONE type elongatedcoil real(kind=db):: axialdims(2) real(kind=db):: radialdims(2) real(kind=db):: center(2) integer:: nbsubcoils(2) real(kind=db):: current real(kind=db):: Nturns end type elongatedcoil type(elongatedcoil), allocatable:: the_coils(:) contains !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Initialize the magnetic field from the standard input !> @param[in] fileid file identifier of the input corresponding to the run parameters !> @param[in] cstep current simulation step number !--------------------------------------------------------------------------- subroutine magnet_init(fileid,cstep) use basic, ONLY:B0,rnorm,bnorm,rgrid,zgrid,width,Rcurv,mpirank,Athet,Br,Bz,bscaling,magnetfile use mpi IMPLICIT NONE integer:: fileid, cstep, istat, ierr integer :: magfiletype = 0 character(256):: line='' NAMELIST /magnetparams/ magnetfile, magfiletype, bscaling ! read the input parameters from file Rewind(fileid) READ(fileid, magnetparams, iostat=istat) ! save the parameters on output IF(mpirank .eq. 0) THEN WRITE(*, magnetparams) END IF if (istat.gt.0) then backspace(fileid) read(fileid,fmt='(A)') line write(*,'(A)') & 'Invalid line in magnetparams: '//trim(line) call MPI_Abort(MPI_COMM_WORLD, -1, ierr) stop end if if(magfiletype .eq. 1)then call read_magfile_txt(magnetfile,the_coils) call compute_B_on_grid(Athet,Br,Bz,rgrid*rnorm,zgrid*rnorm) else if(magfiletype .eq. 0 .and. len_trim(magnetfile) .lt. 1)then call magmirror(Athet,Br,Bz,rgrid,zgrid,rnorm,B0,Rcurv,width) else CALL load_mag_from_h5(magnetfile,Athet,Br,Bz,rgrid,zgrid, B0, rnorm,bscaling) end if ! We normalise the magnetic fields Br=Br/bnorm Bz=Bz/bnorm end subroutine !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Loads the magnetic field defined in the .h5 file at location magfile !> @param[in] magfile filname of .h5 file containing the definitions of A and B !--------------------------------------------------------------------------- SUBROUTINE load_mag_from_h5(magfile, Athet,Br,Bz,rgrid,zgrid, B0, rnorm,bscaling) USE constants, ONLY: Pi USE futils USE bsplines CHARACTER(LEN=*), INTENT(IN):: magfile REAL(kind=db), ALLOCATABLE :: magr(:), magz(:) REAL(kind=db), ALLOCATABLE :: tempBr(:, :), tempBz(:, :), tempAthet(:, :) real(kind=db), allocatable:: c(:,:) real(kind=db)::B0,rnorm Integer:: bscaling, nr, nz,i real(kind=db):: Athet(:), Br(:), Bz(:) Real(kind=db):: zgrid(0:), rgrid(0:) Real(kind=db), ALLOCATABLE:: vec1(:), vec2(:) type(spline2d):: Maginterpolation REAL(kind=db) :: maxB INTEGER :: magfid, dims(2) LOGICAL:: B_is_saved INTEGER :: magn(2), magrank CALL openf(trim(magfile), magfid, 'r', real_prec='d') CALL getdims(magfid, '/mag/Athet', magrank, magn) nr=size(rgrid,1)-1 nz=size(zgrid,1)-1 ! Auxiliary vectors ALLOCATE(vec1((nz+1)*(nr+1)),vec2((nr+1)*(nz+1))) DO i=0,nr vec1(i*(nz+1)+1:(i+1)*(nz+1))=zgrid!(0:nz) vec2(i*(nz+1)+1:(i+1)*(nz+1))=rgrid(i) END DO ALLOCATE (magr(magn(2)), magz(magn(1))) ALLOCATE (tempAthet(magn(1), magn(2)), tempBr(magn(1), magn(2)), tempBz(magn(1), magn(2))) ! Read r and z coordinates for the definition of A_\thet, and B CALL getarr(magfid, '/mag/r', magr) CALL getarr(magfid, '/mag/z', magz) CALL getarr(magfid, '/mag/Athet', tempAthet) IF (isdataset(magfid, '/mag/Br') .and. isdataset(magfid, '/mag/Bz')) THEN CALL getarr(magfid, '/mag/Br', tempBr) CALL getarr(magfid, '/mag/Bz', tempBz) IF(bscaling .gt. 0) then maxB=sqrt(maxval(tempBr**2+tempBz**2)) tempBr=tempBr/maxB*B0 tempBz=tempBz/maxB*B0 end if B_is_saved = .true. ELSE B_is_saved = .false. END IF magz=magz/rnorm magr=magr/rnorm CALL set_splcoef((/3,3/),magz,magr,Maginterpolation) call get_dim(Maginterpolation,dims) ! Interpolation of the magnetic potential vector allocate(c(dims(1),dims(2))) call get_splcoef(Maginterpolation,tempAthet, c) CALL gridval(Maginterpolation,vec1,vec2, Athet ,(/0,0/),c) - if(B_is_saved == .true.)then + if(B_is_saved .eqv. .true.)then ! Interpolation of the Axial magnetic field call get_splcoef(Maginterpolation,tempBz, c) CALL gridval(Maginterpolation,vec1,vec2, Bz ,(/0,0/),c) ! Interpolation of the radial magnetic field call get_splcoef(Maginterpolation,tempBr, c) CALL gridval(Maginterpolation,vec1,vec2, Br ,(/0,0/),c) else CALL gridval(Maginterpolation,vec1,vec2, Br,(/1,0/)) Br=-Br CALL gridval(Maginterpolation,vec1,vec2, Bz,(/0,1/)) Bz=Bz+Athet/vec2 end if if( bscaling .lt. 0 ) then maxB = maxval(sqrt(Bz**2 + Br**2)) Bz = Bz/maxB*B0 Br = Br/maxB*B0 end if CALL closef(magfid) deallocate(c) call destroy_SP(Maginterpolation) END SUBROUTINE load_mag_from_h5 !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Computes the magnetic field on the grid according to a magnetic mirror !> @param[out] Athet magnetic potential vector !> @param[out] Br radial magnetic field on the grid !> @param[out] Bz Axial magnetic field on the grid !> @param[in] rgrid radial grid points !> @param[in] zgrid axial grid points !> @param[in] rnorm normalisation distance constant !> @param[in] B0 maximum mag field amplitude on axis !> @param[in] Rcurv magnetic mirror R=max(B)/min(B) parameter !> @param[in] width distance between the two mirror coils [m] !> !--------------------------------------------------------------------------- subroutine magmirror(Athet,Br,Bz,rgrid,zgrid,rnorm,B0,Rcurv,width) USE constants, ONLY: Pi real(kind=db):: Athet(:), Br(:), Bz(:) REAL(kind=db) :: rg, zg, halfLz, rnorm, Rcurv, MirrorRatio, width, B0 Real(kind=db):: zgrid(:), rgrid(:) INTEGER :: i, rindex, nr, nz nr=size(rgrid,1)-1 nz=size(zgrid,1)-1 halfLz = (zgrid(nz) + zgrid(0))/2 MirrorRatio = (Rcurv - 1)/(Rcurv + 1) DO i = 1, (nr + 1)*(nz + 1) rindex = (i - 1)/(nz + 1) rg = rgrid(rindex) zg = zgrid(i - rindex*(nz + 1) - 1) - halfLz Br(i) = -B0*MirrorRatio*SIN(2*pi*zg/width*rnorm)*bessi1(2*pi*rg/width*rnorm) Bz(i) = B0*(1 - MirrorRatio*COS(2*pi*zg/width*rnorm)*bessi0(2*pi*rg/width*rnorm)) Athet(i) = 0.5*B0*(rg*rnorm - width/pi*MirrorRatio*bessi1(2*pi*rg/width*rnorm)*COS(2*pi*zg/width*rnorm)) END DO end subroutine !________________________________________________________________________________ !Modified Bessel functions of the first kind of the zero order FUNCTION bessi0(x) REAL(kind=db) :: bessi0, x REAL(kind=db) :: ax REAL(kind=db) p1, p2, p3, p4, p5, p6, p7, q1, q2, q3, q4, q5, q6, q7, q8, q9, y SAVE p1, p2, p3, p4, p5, p6, p7, q1, q2, q3, q4, q5, q6, q7, q8, q9 DATA p1, p2, p3, p4, p5, p6, p7/1.0d0, 3.5156229d0, 3.0899424d0, 1.2067492d0, 0.2659732d0, 0.360768d-1, 0.45813d-2/ DATA q1, q2, q3, q4, q5, q6, q7, q8, q9/0.39894228d0, 0.1328592d-1, 0.225319d-2, -0.157565d-2, 0.916281d-2, & & -0.2057706d-1, 0.2635537d-1, -0.1647633d-1, 0.392377d-2/ if (abs(x) .lt. 3.75) then y = (x/3.75)**2 bessi0 = p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7))))) else ax = abs(x) y = 3.75/ax bessi0 = (exp(ax)/sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9)))))))) end if return END FUNCTION bessi0 !________________________________________________________________________________ !Modified Bessel functions of the first kind of the first order FUNCTION bessi1(x) REAL(kind=db) :: bessi1, x REAL(kind=db) :: ax REAL(kind=db) p1, p2, p3, p4, p5, p6, p7, q1, q2, q3, q4, q5, q6, q7, q8, q9, y SAVE p1, p2, p3, p4, p5, p6, p7, q1, q2, q3, q4, q5, q6, q7, q8, q9 DATA p1, p2, p3, p4, p5, p6, p7/0.5d0, 0.87890594d0, 0.51498869d0, 0.15084934d0, 0.2658733d-1, 0.301532d-2, 0.32411d-3/ DATA q1, q2, q3, q4, q5, q6, q7, q8, q9/0.39894228d0, -0.3988024d-1, -0.362018d-2, 0.163801d-2, -0.1031555d-1, & & 0.2282967d-1, -0.2895312d-1, 0.1787654d-1, -0.420059d-2/ if (abs(x) .lt. 3.75D0) then y = (x/3.75D0)**2 bessi1 = x*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7)))))) else ax = abs(x) y = 3.75D0/ax bessi1 = (exp(ax)/sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9)))))))) if (x .lt. 0.) bessi1 = -bessi1 end if return END FUNCTION bessi1 ! calculation of the magnetic potential vector and magnetic fields ! using elliptic integrals ! see Smythe "static and dynamic electricity" McGraw-Hill, third edition p.290 subroutine compute_B_on_grid(A,Br,Bz,r,z) use elliptic use constants implicit none real(kind=db):: A(:), Br(:), Bz(:) real(kind=db):: r(:), z(:) integer:: i,j,k, l, n,linindex, nr,nz, ncoils Real(kind=db):: kappa2, rho, distz, dr, dz Real(kind=db):: r_wire, z_wire Real(kind=db):: r_cen, z_cen, u, v, w, t Real(kind=db):: EE, EK Real(kind=db):: Isubcoil, coilsurf, Itot, currfac nr=size(r,1) nz=size(z,1) ncoils=size(the_coils,1) !$OMP PARALLEL DO PRIVATE(j,i,linindex,k,dr,dz,z_cen,r_cen,coilsurf,Itot,Isubcoil,l,z_wire,n,r_wire,kappa2,EE,EK, u,v,w,t,distz, currfac) Do j=1,nr Do i=1,nz linindex=i+(j-1)*size(z,1) A(linindex)=0.0_db Bz(linindex)=0.0_db Br(linindex)=0.0_db Do k=1,ncoils dr = (the_coils(k)%radialdims(2)-the_coils(k)%radialdims(1)) dz = (the_coils(k)%axialdims(2)-the_coils(k)%axialdims(1)) r_cen = (the_coils(k)%radialdims(2)+the_coils(k)%radialdims(1))/2 z_cen = (the_coils(k)%axialdims(2)+the_coils(k)%axialdims(1))/2 coilsurf = dr * dz; Itot = the_coils(k)%current*the_coils(k)%Nturns Isubcoil = Itot/(the_coils(k)%nbsubcoils(1)*the_coils(k)%nbsubcoils(2)) currfac = mu_0*Isubcoil/pi Do l=1,the_coils(k)%nbsubcoils(1) z_wire = (z_cen - dz/2) + (dz/the_coils(k)%nbsubcoils(1))*(l-0.5) distz = (z(i)-z_wire) Do n=1,the_coils(k)%nbsubcoils(2) r_wire=(r_cen - dr/2) + (dr/the_coils(k)%nbsubcoils(2))*(n-0.5); u=(r_wire+r(j))**2+distz**2 v=(r_wire-r(j))**2+distz**2 w=r_wire**2+r(j)**2+distz**2 t=r_wire**2-r(j)**2-distz**2 kappa2=4*r_wire*r(j)/u EK=elliptic_fk(sqrt(kappa2)) ! complete Elliptic integral of the first kind EE=elliptic_ek(sqrt(kappa2)) ! complete Elliptic integral of the second kind if(r(j).gt.0.0_db)then A(linindex)=A(linindex)+currfac*sqrt(r_wire/(r(j)*kappa2))*& & ((1-kappa2/2)*EK-EE) Br(linindex)=Br(linindex)+currfac/2* & & distz/(r(j)*sqrt(u))*(-EK+w/v*EE) end if Bz(linindex)=Bz(linindex)+currfac/2* & & 1/sqrt(u)*(EK+t/v*EE) end do end do end do end do !WRITE(*,*)"id done",j,'over',nr end do !$OMP END PARALLEL DO end subroutine subroutine read_magfile_txt(magnetfile,the_coils) type(elongatedcoil), allocatable:: the_coils(:) character(128):: magnetfile INTEGER:: nb_coils INTEGER :: lu_magfile=999 INTEGER:: i, openerr, reason CHARACTER(len=256) :: header real(kind=db):: z1,z2, r1, r2, current integer:: Nr, Nz real(kind=db):: Nturns if(allocated(the_coils)) deallocate(the_coils) nb_coils=0 OPEN(UNIT=lu_magfile,FILE=trim(magnetfile),ACTION='READ',IOSTAT=openerr) header=' ' IF(openerr .ne. 0) THEN CLOSE(unit=lu_magfile) RETURN END IF ! The magnet table is defined as a 8 column (zmin zmax rmin rmax nturns Nz Nr current) DO WHILE(.true.) READ(lu_magfile,'(a)',IOSTAT=reason) header header=adjustl(header) if(reason .lt. 0 ) exit ! We reached end of file if( header(1:1) .ne. '!'.and. len_trim(header) .gt. 1) then READ(header,*,IOSTAT=reason) z1,z2, r1, r2, nturns, Nz, Nr, current if(reason .lt. 0 ) then WRITE(*,*) "Error in magnet definition of file: ",magnetfile,"\n with input ", header exit end if nb_coils=nb_coils+1 end if END DO allocate(the_coils(nb_coils)) REWIND(lu_magfile) ! The magnet table is defined as a 8 column (zmin zmax rmin rmax nturns Nz Nr current) i=1 write(*,*) "number of found coils", nb_coils DO WHILE(i .le. nb_coils) READ(lu_magfile,'(a)',IOSTAT=reason) header header=adjustl(header) !Write(*,*) i,' ',trim(header) if(reason .lt. 0 ) exit ! We reached end of file if( header(1:1) .ne. '!' .and. len_trim(header) .gt. 1) then READ(header,*,IOSTAT=reason) z1,z2, r1, r2, nturns, Nz, Nr, current the_coils(i)%axialdims=(/z1,z2/) the_coils(i)%radialdims=(/r1,r2/) the_coils(i)%Nturns=nturns the_coils(i)%nbsubcoils=(/Nz,Nr/) the_coils(i)%current=current !Write(*,*) i,' ',the_coils(i)%axialdims, the_coils(i)%radialdims, the_coils(i)%Nturns,the_coils(i)%nbsubcoils,the_coils(i)%current i=i+1 end if END DO CLOSE(unit=lu_magfile) end subroutine -end module magnet \ No newline at end of file +end module magnet diff --git a/src/neutcol_mod.f90 b/src/neutcol_mod.f90 index 91979be..9016868 100644 --- a/src/neutcol_mod.f90 +++ b/src/neutcol_mod.f90 @@ -1,574 +1,574 @@ !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: neutcol ! !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> Module responsible for handling the electron-neutral collisions and creating electrons !> by ionisation Based on the paper by Birdsall 1991 and Sengupta et al. !------------------------------------------------------------------------------ module neutcol USE constants IMPLICIT NONE private LOGICAL, SAVE :: nlcol=.false. !< Flag to activate or not electron neutral collisions LOGICAL :: nlmaxwellio=.false. !< Flag to define how ionised electrons are created (physically or according to maxwellian) INTEGER :: itcol = 1 !< number of dt between each evaluation of neutcol_step Real(kind=db) :: neutdens=2.4e16 !< Neutral particle density in m-3 Real(kind=db) :: neuttemp=300 !< Neutral particle temperature in K Real(kind=db) :: neutpressure !< Neutral particle pressure in mbar Real(kind=db) :: scatter_fac = 24.2 !< Energy scattering factor for the considered gas (here for Ne) [eV] see Opal 1971 https://doi.org/10.1063/1.1676707 real(kind=db) :: Eion = 21.56 !< Ionisation energy (eV) (here for Ne) Real(kind=db) :: E0 = 27.21 !< Atomic unit of energy used for calculation of deviation angles [eV] Real(kind=db) :: dE = 1 !< resolution for the computation of the cross-sections Real(kind=db) :: Emax = 1 !< resolution for the computation of the cross-sections real(kind=db) :: collfactor !< Normalised collision factor (n_n \delta t) INTEGER :: nb_io_cross=0 Real(kind=db), ALLOCATABLE :: io_cross_sec(:,:) !< Ionisation cross-section table Real(kind=db), ALLOCATABLE :: io_cross_sec_lin(:,:) !< Ionisation cross-section table for linear interpolation Real(kind=db), ALLOCATABLE :: io_growth_cross_sec(:) !< Ionisation exponential fitting factor INTEGER :: nb_ela_cross=0 Real(kind=db), ALLOCATABLE :: ela_cross_sec(:,:) !< Elastic collision cross section table Real(kind=db), ALLOCATABLE :: ela_cross_sec_lin(:,:) !< Elastic collision cross section table for linear interpolation Real(kind=db), ALLOCATABLE :: ela_growth_cross_sec(:) !< Elastic collision exponential fitting factor Real(kind=db) :: Escale !< Energy normalisation factor used to reduce computation costs CHARACTER(len=128) :: io_cross_sec_file='' CHARACTER(len=128) :: ela_cross_sec_file='' Real(kind=db) :: etemp=22000 !< In case of nlmaxwelio, defines the temperature of created electrons [K] Real(kind=db) :: vth !< In case of nlmaxwelio, defines the normalised thermal velocity of created electrons LOGICAL :: nldragio=.true. !< Set if inpinging electrons are affected by ionising collisions INTEGER :: species(2) !< species(1) contains the specie index in plist which stores the colliding particles, species(2) stores the specie index for the released ion. LOGICAL :: isotropic = .false. !< is the scattering angle isotropic type coll_process Real(kind=db):: Emin !< minimum energy in EV for collision to happen Real(kind=db):: Emax !< maximum energy in EV for linear interpolation Real(kind=db):: deltaE !< energy step in ev for interpolation Real(kind=db), allocatable:: cross_sec(:) !< cross-section as read from file Real(kind=db), allocatable:: E_levels(:) !< energy values as read from file Real(kind=db), allocatable:: growth_parameter(:) !< collision exponential fitting factor Real(kind=db), allocatable:: lin_E(:) !< linear energy levels for linear interpolation Real(kind=db), allocatable:: lin_cross(:) !< linear cross-sections for linear interpolation end type NAMELIST /neutcolparams/ neutdens, Eion, & & scatter_fac, nlcol, io_cross_sec_file, ela_cross_sec_file, nlmaxwellio, etemp, & & nldragio, itcol, species, isotropic PUBLIC:: neutcol_init, neutcol_step, neutcol_diag, itcol, neutdens PROCEDURE(rotate_vel), POINTER:: change_dir => NULL()!< Function evaluating the weight for Dirichelt boundary conditions ABSTRACT INTERFACE SUBROUTINE rotate_vel(Ur, Uthet, Uz, coschi, thet) use constants real(kind=db), INTENT(INOUT):: Ur, uthet, uz, coschi, thet END SUBROUTINE end interface CONTAINS subroutine neutcol_init(lu_in, p) use mpi Use basic, only: mpirank, dt, nlclassical,rnorm, vnorm Use beam, only: particles Use constants implicit none INTEGER, INTENT(IN) :: lu_in TYPE(particles) :: p INTEGER:: ierr, istat, i character(len=1000) :: line real(kind=db):: xsi species(1)=1 species(2)=-1 Rewind(lu_in) READ(lu_in, neutcolparams, iostat=istat) if (istat.gt.0) then backspace(lu_in) read(lu_in,fmt='(A)') line write(*,'(A)') & 'Invalid line in neutcolparams: '//trim(line) call MPI_Abort(MPI_COMM_WORLD, -1, ierr) stop end if IF(mpirank .eq. 0) THEN WRITE(*, neutcolparams) END IF if(.not. nlcol) return if(nlclassical)THEN Escale=0.5*p%m/elchar*vlight**2 else Escale=p%m*vlight**2/elchar end if if (nlmaxwellio) vth=sqrt(kb*etemp/p%m)/vnorm if(io_cross_sec_file .ne.'') then call read_cross_sec(io_cross_sec_file,io_cross_sec, nb_io_cross) if(nb_io_cross .gt. 0) then allocate(io_growth_cross_sec(nb_io_cross-1)) ! Normalisations io_cross_sec(:,2)=io_cross_sec(:,2)/rnorm**2 ! Precomputing of exponential fitting factor for faster execution io_growth_cross_sec=log(io_cross_sec(2:nb_io_cross,2)/io_cross_sec(1:nb_io_cross-1,2))/ & & log(io_cross_sec(2:nb_io_cross,1)/io_cross_sec(1:nb_io_cross-1,1)) end if end if if(ela_cross_sec_file .ne.'') then call read_cross_sec(ela_cross_sec_file,ela_cross_sec, nb_ela_cross) if(nb_ela_cross .gt. 0) then allocate(ela_growth_cross_sec(nb_ela_cross-1)) ! Normalisations ela_cross_sec(:,2)=ela_cross_sec(:,2)/rnorm**2 if(.not. isotropic) then do i=1,nb_ela_cross xsi=ela_cross_sec(i,1)/(0.25*E0+ela_cross_sec(i,1)) ela_cross_sec(i,2)=ela_cross_sec(i,2)*(2*xsi**2)/((1-xsi)*((1+xsi)*log((1+xsi)/(1-xsi))-2*xsi)) end do end if ! Precomputing of exponential fitting factor for faster execution ela_growth_cross_sec=log(ela_cross_sec(2:nb_ela_cross,2)/ela_cross_sec(1:nb_ela_cross-1,2))/ & & log(ela_cross_sec(2:nb_ela_cross,1)/ela_cross_sec(1:nb_ela_cross-1,1)) end if END IF nlcol=nlcol .and. (allocated(io_cross_sec) .or. allocated(ela_cross_sec)) ! Collision factor depending on neutral gas parameters collfactor=neutdens*dt*rnorm**3*itcol neutpressure=neutdens*kb*300/100 if (.not. isotropic)then change_dir=> rotate else change_dir=> scatter end if end subroutine neutcol_init Subroutine neutcol_diag(File_handle, str, vnorm) use mpi Use futils Integer:: File_handle Real(kind=db):: vnorm Character(len=*):: str CHARACTER(len=256):: grpname Integer:: ierr, mpirank,i Real(kind=db)::xsi Real(kind=db),allocatable:: nonisotropic_rescale(:,:) CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) IF(mpirank .eq. 0 .and. nlcol) THEN Write(grpname,'(a,a)') trim(str),"/neutcol" If(.not. isgroup(File_handle, trim(grpname))) THEN CALL creatg(File_handle, trim(grpname)) END IF Call attach(File_handle, trim(grpname), "neutdens", neutdens) Call attach(File_handle, trim(grpname), "neuttemp", neuttemp) Call attach(File_handle, trim(grpname), "neutpressure", neutpressure) Call attach(File_handle, trim(grpname), "scatter_fac", scatter_fac) Call attach(File_handle, trim(grpname), "Eion", Eion) Call attach(File_handle, trim(grpname), "E0", E0) Call attach(File_handle, trim(grpname), "Escale", Escale) Call putarr(File_handle,trim(grpname)//"species", species) if (allocated(io_cross_sec)) Call putarr(File_handle, trim(grpname)//"/io_cross_sec", io_cross_sec) allocate(nonisotropic_rescale(nb_ela_cross,2)) nonisotropic_rescale=1 if(.not. isotropic) then do i=1,nb_ela_cross xsi=ela_cross_sec(i,1)/(0.25*E0+ela_cross_sec(i,1)) nonisotropic_rescale(i,2)=(2*xsi**2)/((1-xsi)*((1+xsi)*log((1+xsi)/(1-xsi))-2*xsi)) end do end if if (allocated(ela_cross_sec)) Call putarr(File_handle, trim(grpname)//"/ela_cross_sec", ela_cross_sec/nonisotropic_rescale) END IF End subroutine neutcol_diag !------------------------------------------------------------- !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Simulates the elastic and ionising collisions for each particles in plist(species(1)) ! !> @param [inout] plist list of particle species considered in the code !--------------------------------------------------------------------------- SUBROUTINE neutcol_step(plist) ! - USE random + USE random_mod USE beam USE omp_lib USE basic, ONLY: nlclassical USE distrib, ONLY: lodgaus type(particles), TARGET::plist(:) type(particles),pointer::p INTEGER:: i, omp_thread, num_threads, j, nbcolls_ela, nbcolls_io real(kind=db):: Rand(5) real(kind=db):: v2, v, ek, Everif, es, cosChi, thet, sig_io, sig_ela, vfact, xsi type(linked_part_row):: ins_p type(linked_part), POINTER:: created real(kind=db):: collisionfact,nucol(3),vinit(3),vend(3) p=>plist(species(1)) if(.not. nlcol .or. p%nploc .le. 0) return num_threads=omp_get_max_threads() nbcolls_ela=0 nbcolls_io=0 nucol=0 !!$OMP private(collisionfact,i,omp_thread,Rand,v2,ek,sig_io,sig_ela,es,coschi,thet,vfact, created, v, everif,xsi,vinit,vend)!, reduction(+:nbcolls_ela,nbcolls_io, nucol) omp_thread=omp_get_thread_num()+1 !omp_thread=1 allocate(ins_p%start) ins_p%n=0 created=>ins_p%start !$OMP DO DO i=1,p%Nploc !for each particle CALL random_array(Rand,1,ran_index(omp_thread),ran_array(:,omp_thread)) ! we calculate the kinetic energy and norm of the velocity v2=(p%U(1,i)**2+p%U(2,i)**2+p%U(3,i)**2) if(nlclassical) THEN ek=v2*escale v=sqrt(v2) vinit=p%U(:,i) ! (/p%UR(i),p%UTHET(i),p%UZ(i)/) ELSE ek=(p%gamma(i)-1)*escale v=sqrt(v2)/p%gamma(i) vinit=p%U(:,i)/p%gamma(i)!(/p%UR(i),p%UTHET(i),p%UZ(i)/)/p%gamma(i) end if sig_io=0 sig_ela=0 ! computes the ionisation and elastic collision cross-sections at this kinetic energy ! The ionisation event can only occur if the incoming electron energy is above the binding energy if (ek .gt. Eion .and. nb_io_cross .gt. 1) then sig_io=sig_fit(io_cross_sec,io_growth_cross_sec,ek, nb_io_cross) end if if (nb_ela_cross .gt. 1) then sig_ela=sig_fit(ela_cross_sec,ela_growth_cross_sec,ek, nb_ela_cross) end if collisionfact=1-exp(-collfactor*(sig_io+sig_ela)*v) ! If we have a collision event if (Rand(1) .lt.collisionfact) THEN CALL random_array(Rand,1,ran_index(omp_thread),ran_array(:,omp_thread)) ! Check if elastic or ionising event is happening IF(Rand(1).gt. sig_ela/(sig_io+sig_ela)) THEN ! An ionisation collision happened and we create the necessary electron ! prepare the memory for the released electron ins_p%n=ins_p%n+1 allocate(created%next) created%next%prev=>created ! Fill created particle new position created%p%pos=p%pos(:,i)!(/p%R(i), p%THET(i), p%Z(i)/) IF( nlmaxwellio ) THEN ! the new electron velocity is defined according to a Maxwellian CALL lodgaus(0, Rand(1:3)) ! get random velocity created%p%U=vth*Rand(1:3) ELSE CALL random_array(Rand,3,ran_index(omp_thread),ran_array(:,omp_thread)) ! Compute created electron energy Es=scatter_fac*tan(Rand(1)*atan((Ek-Eion)/(2*scatter_fac))) ! Compute scattering angles for created electron if (isotropic) then coschi=cos(Rand(2)*pi) else cosChi=1-2*Rand(2)/(1+8*Es/E0*(1-Rand(2))) end if thet=Rand(3)*2*pi if(nlclassical)THEN ! new velocity factor for created particle vfact=sqrt(Es/Ek) ELSE ! new velocity factor for created particle vfact=sqrt(Es*(Es+2*Escale)/(Ek*(Ek+2*Escale))) END IF ! Fill created particle velocity created%p%U=vfact*p%U(:,i)!(/p%UR(i),p%UTHET(i),p%UZ(i)/) ! rotate the velocity vector due to the collision call change_dir(created%p%U(1),created%p%U(2), created%p%U(3), coschi, thet) END IF vend=created%p%U if(nlclassical)THEN ! Lorentz factor for created particle created%p%gamma=1.0 ELSE ! Lorentz factor for created particle created%p%gamma=sqrt(1+created%p%U(1)**2+created%p%U(2)**2+created%p%U(3)**2) vend=vend/created%p%gamma END IF ! We prepare the next created particle ins_p%end=>created created=>created%next ! We keep track of what changed nbcolls_io=nbcolls_io+1 nucol=nucol-vend/vinit ! If we want the incoming electron to be scattered, we need to compute ! its new kinetic energy if (nldragio) THEN ! We store the lossed energy in pot for keeping track of energy conservation created%prev%p%pot=Eion+Es CALL random_array(Rand,2,ran_index(omp_thread),ran_array(:,omp_thread)) Es=Ek-Eion-Es if(nlclassical)THEN ! new velocity factor for scattered particle vfact=sqrt(Es/Ek) ELSE ! new velocity factor for scattered particle vfact=sqrt(Es*(Es+2*Escale)/(Ek*(Ek+2*Escale))) END IF ELSE CYCLE END IF ELSE ! An elastic collision event happens CALL random_array(Rand,2,ran_index(omp_thread),ran_array(:,omp_thread)) Es=Ek vfact=1.0 nbcolls_ela=nbcolls_ela+1 END IF ! We calculate the scattered velocity angle for the scattered electron if (isotropic) then coschi=cos(Rand(1)*pi) else cosChi=1-2*Rand(1)/(1+8*Es/E0*(1-Rand(1))) end if thet=Rand(2)*2*pi ! Change the incident electron velocity direction and amplitude if necessary p%U(:,i)=p%U(:,i)*vfact !p%UTHET(i)=p%UTHET(i)*vfact !p%UZ(i)=p%UZ(i)*vfact call change_dir(p%U(1,i),p%U(2,i), p%U(3,i), coschi, thet) if(nlclassical) THEN vend=p%U(:,i)!(/p%UR(i),p%UTHET(i),p%UZ(i)/) ELSE p%gamma(i)=sqrt(1+p%U(1,i)**2+p%U(2,i)**2+p%U(3,i)**2) vend=p%U(:,i)/p%gamma(i)!(/p%UR(i),p%UTHET(i),p%UZ(i)/)/p%gamma(i) END IF nucol=nucol+1-vend/vinit END IF END DO !$OMP END DO NOWAIT !$OMP BARRIER ! clean up the memory after the loop if(associated(created%prev)) then created=>created%prev ins_p%end=>created deallocate(created%next) else deallocate(ins_p%start) end if !!$OMP END PARALLEL ! We collect all created particules into one linked list for easier insertion in plist !Do i=1,num_threads !$OMP CRITICAL (insertions) if(species(2).gt.0.and.ins_p%n .gt.0) then CALL add_created_part(plist(species(2)), ins_p, .false.,.true.) end if !$OMP END CRITICAL (insertions) !$OMP CRITICAL (insertelectrons) if(ins_p%n .gt.0) then CALL add_created_part(plist(species(1)),ins_p,.true.,.false.) !exit end if !$OMP END CRITICAL (insertelectrons) !end do !$OMP CRITICAL(addcolls) p%nbcolls=p%nbcolls+(/nbcolls_io, nbcolls_ela/) p%nudcol=p%nudcol+nucol !$OMP END CRITICAL(addcolls) !Write(*,*)"mpirank: ", mpirank, " Nb colls ela, io: ",nbcolls_ela, nbcolls_io ! END SUBROUTINE neutcol_step FUNCTION sig_fit(sig_vec,growth_vec,ek,nb_cross) use distrib, ONLY: closest real(kind=db)::sig_fit, ek real(kind=db):: sig_vec(:,:), growth_vec(:) Integer:: k, nb_cross sig_fit=0 k=closest(sig_vec(:,1),ek, nb_cross-1) if(k.lt.1) return !sig_fit=(sig_vec(k,1)-sig_vec(k-1,1))/(sig_vec(k,2)-sig_vec(k-1,2))*(sig_vec(k,2)-ek)+sig_vec(k-1,1) ! Exponential fitting relevant at high energies sig_fit=sig_vec(k,2)*(ek/sig_vec(k,1))**growth_vec(k) END FUNCTION sig_fit FUNCTION sig_fit_lin(sig_vec,growth_vec,ek,nb_cross) use distrib, ONLY: closest real(kind=db)::sig_fit_lin, ek real(kind=db):: sig_vec(:,:), growth_vec(:) Integer:: k, nb_cross sig_fit_lin=0 k=(k-sig_vec(k,1))/dE if(k.lt.1) return !sig_fit=(sig_vec(k,1)-sig_vec(k-1,1))/(sig_vec(k,2)-sig_vec(k-1,2))*(sig_vec(k,2)-ek)+sig_vec(k-1,1) ! Exponential fitting relevant at high energies sig_fit_lin=sig_vec(k,2)+(ek-sig_vec(k,1))/(sig_vec(k+1,1)-sig_vec(k,1))*(sig_vec(k+1,2)-sig_vec(k,2)) END FUNCTION sig_fit_lin SUBROUTINE rotate(Ur, Uthet, Uz, coschi, thet) real(kind=db), INTENT(INOUT):: Ur, uthet, uz, coschi, thet real(kind=db):: norm, perp(3), U(3), U0(3) real(kind=db):: sinchi, sinthet, costhet Integer :: iperp1,iperp2 U0=(/Ur,Uthet,Uz/) norm=sqrt(sum(U0**2)) U=U0/norm ! Find a vector perpendicular to U for chi rotation ! find the direction with maximum amplitude perp=(/1,1,1/) iperp1=maxloc(abs(U),1) ! find second direction with next max amplitude perp(iperp1)=0 iperp2=maxloc(abs(perp*U),1) perp=0 perp(iperp2)=U(iperp1) perp(iperp1)=-U(iperp2) ! Normalise the rotation vector perp=perp/sqrt(sum(perp**2)) ! Compute sinus and cosinus for rotation sinchi=sqrt(1-coschi**2) costhet=cos(thet) sinthet=sin(thet) ! Rotation of angle chi around perp Ur = (coschi+perp(1)**2*(1-coschi))*U0(1) + (perp(1)*perp(2)*(1-coschi)-perp(3)*sinchi)*U0(2) + (perp(1)*perp(3)*(1-coschi) + perp(2)*sinchi)*U0(3) Uthet = (perp(1)*perp(2)*(1-coschi)+perp(3)*sinchi)*U0(1) + (coschi + perp(2)**2*(1-coschi))*U0(2) +(perp(2)*perp(3)*(1-coschi)-perp(1)*sinchi)*U0(3) Uz = (perp(1)*perp(3)*(1-coschi)-perp(2)*sinchi)*U0(1) +(perp(3)*perp(2)*(1-coschi)+perp(1)*sinchi)*U0(2) +( coschi+perp(3)**2*(1-coschi))*U0(3) U0 =(/Ur,Uthet,Uz/) ! second rotation according to uniform distribution ! Rotation of angle theta around U Ur = (costhet+U(1)**2*(1-costhet))*U0(1) + (U(1)*U(2)*(1-costhet) - U(3)*sinthet)*U0(2) + (U(1)*U(3)*(1-costhet)+U(2)*sinthet)*U0(3) Uthet = (U(2)*U(1)*(1-costhet)+U(3)*sinthet)*U0(1) + (costhet + U(2)**2*(1-costhet))*U0(2) + (U(2)*U(3)*(1-costhet)-U(1)*sinthet)*U0(3) Uz = (U(3)*U(1)*(1-costhet) - U(2)*sinthet)*U0(1) + (U(3)*U(2)*(1-costhet)+U(1)*sinthet)*U0(2) + (costhet +U(3)**2*(1-costhet))*U0(3) !normf=sqrt(Ur**2+Uthet**2+Uz**2) !if(abs(norm-normf)/norm .gt. 1e-14) WRITE(*,*) "Error in rotate the norm of v changed" END SUBROUTINE rotate SUBROUTINE scatter(Ur, Uthet, Uz, coschi, thet) real(kind=db), INTENT(INOUT):: Ur, uthet, uz, coschi, thet real(kind=db):: norm real(kind=db):: sinchi, sinthet, costhet norm=sqrt(Ur**2+Uz**2+Uthet**2) ! Compute sinus and cosinus for rotation sinchi=sqrt(1-coschi**2) costhet=cos(thet) sinthet=sin(thet) Ur=norm*sinchi*costhet Uthet=norm*sinchi*sinthet Uz=norm*coschi END SUBROUTINE scatter SUBROUTINE read_cross_sec(filename,cross_sec, nb_cross) CHARACTER(len=*) ::filename Real(kind=db), ALLOCATABLE :: cross_sec(:,:) INTEGER:: nb_cross INTEGER :: lu_cross_sec=9999 INTEGER:: i, openerr, reason CHARACTER(len=256) :: header real(kind=db):: t1,t2 nb_cross=0 OPEN(UNIT=lu_cross_sec,FILE=trim(filename),ACTION='READ',IOSTAT=openerr) header=' ' IF(openerr .ne. 0) THEN CLOSE(unit=lu_cross_sec) RETURN END IF ! The cross section table is defined as a two column energy and cross_section DO WHILE(.true.) READ(lu_cross_sec,'(a)',IOSTAT=reason) header header=adjustl(header) if(reason .lt. 0 ) exit ! We reached end of file if( header(1:1) .ne. '!') then READ(header,*) t1, t2 if(t1 .ne. 0 .and. t2.ne. 0) nb_cross=nb_cross+1 end if END DO if (allocated(cross_sec)) deallocate(cross_sec) allocate(cross_sec(nb_cross,2)) REWIND(lu_cross_sec) ! The cross section table is defined as a two column energy and cross_section i=1 DO WHILE(i .le. nb_cross) READ(lu_cross_sec,'(a)',IOSTAT=reason) header header=adjustl(header) if(reason .lt. 0 ) exit ! We reached end of file if( header(1:1) .ne. '!') then READ(header,*) cross_sec(i,1), cross_sec(i,2) if(cross_sec(i,1) .ne. 0 .and. cross_sec(i,2).ne. 0) i=i+1 end if END DO CLOSE(unit=lu_cross_sec) END subroutine read_cross_sec subroutine init_cross_sec_type(cross_sec_data, nb, cross_sec_array, exp_growth_array, Emin, Emax, deltaE) type(coll_process):: cross_sec_data INTEGER:: nb real(kind=db):: cross_sec_array(:,:), exp_growth_array(:) real(kind=db):: Emin, Emax, deltaE allocate(cross_sec_data%cross_sec(nb),cross_sec_data%E_levels(nb)) cross_sec_data%cross_sec=cross_sec_array(:,2) cross_sec_data%E_levels=cross_sec_array(:,1) allocate(cross_sec_data%growth_parameter(nb-1)) cross_sec_data%growth_parameter=exp_growth_array cross_sec_data%Emin=Emin cross_sec_data%Emax=Emax end subroutine end module neutcol