diff --git a/platform/build/make.inc.TUMBLEWEED b/platform/build/make.inc.TUMBLEWEED index c502f5b36..19712de78 100644 --- a/platform/build/make.inc.TUMBLEWEED +++ b/platform/build/make.inc.TUMBLEWEED @@ -14,14 +14,14 @@ NUMAS_PER_NODE=1 # Compilers -FC = mpifort -DDIMATCOPY -std=gnu -fallow-argument-mismatch -fall-intrinsics -fimplicit-none -J $(GACODE_ROOT)/modules #-fPIC +FC = mpifort -DDIMATCOPY -fallow-argument-mismatch -fall-intrinsics -fimplicit-none -J $(GACODE_ROOT)/modules #-fPIC #FC = mpifort -std=f2008 -fall-intrinsics -fimplicit-none -J $(GACODE_ROOT)/modules #-fPIC F77 = mpifort -fimplicit-none -fallow-argument-mismatch FOMP =-fopenmp FMATH =-fdefault-real-8 -fdefault-double-8 -I$(FFTW_INC) -FOPT =-O3 -fexternal-blas -march=native -ffree-line-length-0 +FOPT =-O3 -std=gnu -fexternal-blas -march=native -ffree-line-length-0 # -floop-nest-optimize -floop-parallelize-all -flto -FDEBUG =-Wall -W -fcheck=all -g -fbacktrace -ffpe-trap=invalid,zero,overflow -finit-real=snan # -std=f2003 -fall-intrinsics +FDEBUG =-Wall -W -fcheck=all -Wno-compare-reals -Wno-integer-division -Wno-unused-variable -g -fbacktrace -ffpe-trap=invalid,zero,overflow -finit-real=snan -std=f2008 -fall-intrinsics F2PY = f2py # System math libraries diff --git a/shared/landau/landau.F90 b/shared/landau/landau.F90 index 6532f9b88..a305eed95 100644 --- a/shared/landau/landau.F90 +++ b/shared/landau/landau.F90 @@ -1,26 +1,28 @@ ! gfortran -march=native -O3 -fno-stack-arrays -fimplicit-none -fdefault-real-8 landau.F90 -c !intel: !ifort -stand f15 -warn all -march=native -O3 -heap-arrays 10 -implicitnone -real-size 64 landau.F90 -c -module landau #ifndef __PGI !can't use quad precision with frumpy pgi compiler - real(16), parameter,private :: pi1=atan(1._16)*4 +#define HIREAL real(16) +#define HRS _16 +#else +#define HIREAL real +#define HRS +#endif +module landau + HIREAL, parameter,private :: pi1=atan(1.HRS)*4 ! real, parameter :: intnorm=pi1**(-1.5)*4*pi1 ! int gphys(v) d^3v=1=int g(x) x^2 dx *intnor ! Maybe one should normalise this once for a given xmax and density 1? - real, parameter :: vunit=sqrt(2._16) + real, parameter :: vunit=real(sqrt(2.HRS)) ! v*vunit is v in units of sqrt(T/m) ! v itself is v in units of sqrt(2T/m) ("collision operator units") -#else - real, parameter,private :: pi1=atan(1.)*4 - real, parameter :: vunit=sqrt(2.) -#endif - real, parameter :: intnorm=4/sqrt(pi1) - real, parameter :: normcol=sqrt(8/pi1) + real, parameter :: intnorm=real(4/sqrt(pi1)) + real, parameter :: normcol=real(sqrt(8/pi1)) ! Normalisation of ctest/emat cfield/emat for the self-collisions of the ! reference species without caa or cab. ! Note: must take into account muref^3=sqrt(1/8), since mu^2=m/(2T) - real, parameter :: normcolcgyro=4/sqrt(pi1) + real, parameter :: normcolcgyro=real(4/sqrt(pi1)) integer :: verbose=0 contains @@ -56,7 +58,7 @@ real function tau_ab_inv(mua,mub,nb,za,zb,Ta,Tb) if (present(nb)) nb1=nb if (present(mua)) mua1=mua if (present(mub)) mub1=mub - tau_ab_inv=sqrt(8./pi1)/3*za1**2*zb1**2*nb1/(mua1*Ta1**2)*& + tau_ab_inv=real(sqrt(8./pi1))/3*za1**2*zb1**2*nb1/(mua1*Ta1**2)*& (1+mua1**2*Ta1/(mub1**2*Tb1))/(1+mua1**2/(mub1**2))**1.5 end function tau_ab_inv @@ -71,36 +73,13 @@ subroutine gentestkernel(n,a1,b1,c1,xmax,beta,gp,gw,ngauss,lor_int,dif_int,addcu optional :: t1t2_int ! only calculate that term if it is present. optional :: addcutoff !allocate(lor_int(n,n),dif_int(n,n)) - real p0,p1,pi,q0,q1,qi,r0,r1,ri,s0,s1,si,x,x1,xmax1 + real p0,p1,pi,q0,q1,qi,r0,r1,ri,s0,s1,si,x,xmax1 real f1,f2 integer i,j,k -#ifndef __PGI - real(16) fct1,fct2,fct1lo,fct1hi,fct2lo,fct2hi !inline functions - real(16) w,xmaxx -#else - real fct1,fct2,fct1lo,fct1hi,fct2lo,fct2hi !inline functions - real w,xmaxx -#endif + HIREAL xmaxx,x1 logical t1t2,addcutoff1 real t1,t2 ! for timing -#ifndef __PGI - fct1(w)=.125_16*w**(-3)*(exp(-w*w)*w+(2*w*w-1)*.5_16*sqrt(pi1)*erf(w)) - fct2(w)=.25_16*w**(-1)*(.5_16*sqrt(pi1)*erf(w)-w*exp(-w**2)) - fct2hi(w)=-w**(-1)*((1._16/12)*xmaxx*exp(-xmaxx**2)*(3+2*xmaxx**2)-(sqrt(pi1)/8)*erf(xmaxx)) -#else - fct1(w)=.125*w**(-3)*(exp(-w*w)*w+(2*w*w-1)*.5*sqrt(pi1)*erf(w)) - fct2(w)=.25*w**(-1)*(.5*sqrt(pi1)*erf(w)-w*exp(-w**2)) - fct2hi(w)=-w**(-1)*((1./12)*xmaxx*exp(-xmaxx**2)*(3+2*xmaxx**2)-(sqrt(pi1)/8)*erf(xmaxx)) -#endif - fct1lo(w)=fct1(w)-exp(-xmaxx**2)/6 ! for w1) then do k=1,ngauss x1=xmax*(1+(beta-1)*gp(k)) - x=x1/beta -#ifndef __PGI - f1=fct1hi(real(x1,16)) - f2=fct2hi(real(x1,16)) -#else - f1=fct1hi(x1) - f2=fct2hi(x1) -#endif + x=real(x1)/beta + f1=real(fct1hi(x1)) + f2=real(fct2hi(x1)) do j=1,n if (j==1) then @@ -281,7 +245,7 @@ subroutine gentestkernel(n,a1,b1,c1,xmax,beta,gp,gw,ngauss,lor_int,dif_int,addcu s1=0 si=0 lor_int(j,i)=lor_int(j,i)+pi*ri*f1 - if (t1t2) t1t2_int(j,i)=t1t2_int(j,i)+2*qi*ri*x1*f2 + if (t1t2) t1t2_int(j,i)=t1t2_int(j,i)+2*qi*ri*real(x1)*f2 ! ^^ vgl 1.8.19 (2) do i=2,n s0=si @@ -292,7 +256,7 @@ subroutine gentestkernel(n,a1,b1,c1,xmax,beta,gp,gw,ngauss,lor_int,dif_int,addcu r1=r0 lor_int(j,i)=lor_int(j,i)+pi*ri*f1 dif_int(j,i)=dif_int(j,i)+qi*si*f2 - if (t1t2) t1t2_int(j,i)=t1t2_int(j,i)+2*qi*ri*x1*f2 + if (t1t2) t1t2_int(j,i)=t1t2_int(j,i)+2*qi*ri*real(x1)*f2 enddo enddo enddo @@ -320,8 +284,44 @@ subroutine gentestkernel(n,a1,b1,c1,xmax,beta,gp,gw,ngauss,lor_int,dif_int,addcu ! i.e. we multiply with beta!! if (t1t2) t1t2_int=t1t2_int*beta**(-4) ! ***<--- must be multiplied with (Ta/Tb-1) call cpu_time(t2) - + if (verbose>0) print '(A,G13.3)','gentestkernel took ',t2-t1 + contains + elemental HIREAL function fct1(w) + implicit none + HIREAL,intent(in) :: w + fct1=.125HRS*w**(-3)*(exp(-w*w)*w+(2*w*w-1)*.5HRS*sqrt(pi1)*erf(w)) + end function fct1 + elemental HIREAL function fct2(w) + implicit none + HIREAL,intent(in) :: w + fct2=.25HRS*w**(-1)*(.5HRS*sqrt(pi1)*erf(w)-w*exp(-w**2)) + end function fct2 + elemental HIREAL function fct1lo(w) + implicit none + HIREAL,intent(in) :: w + fct1lo=fct1(w)-exp(-xmaxx**2)/6 ! for w - - w(x)=x**alpha*exp(-sign*x*x) ! inline function for weight - sechb(x)=2.*x/(1.+x*x) - sech(x)=sechb(exp(-x)) - sech2(x)=sech(abs(x))**2 - tanhb(x)=2*x/(1+x) - tanhp1(x)=tanhb(exp(2*x)) ! tanh(x)+1 for x<0 approx. 2*exp(2*x) for x<<-1 - tanhm1(x)=-tanhb(exp(-2*x)) ! tanh(x)-1 for x>0 - if (n==0) return nn=n vb=present(verbose) @@ -182,7 +172,7 @@ subroutine half_hermite_norm(n,xmin,xmax,sign,alpha,a1,b1,c1,a,bsq,logg,verbos& t00=t0 t10=t1 if (vb) then - ngauss=n+xmax**2+xmax*sqrt(-2*log(epsgauss))-log(epsgauss)/3 + ngauss=n+floor(xmax**2+xmax*sqrt(-2*log(epsgauss))-log(epsgauss)/3) ngauss=n+ceiling(10*xmax)+ceiling(alpha/2) print *,'ngauss=',ngauss allocate(gw(ngauss),gp(ngauss),gw2(ngauss*2),gp2(ngauss*2)) @@ -671,6 +661,42 @@ subroutine half_hermite_norm(n,xmin,xmax,sign,alpha,a1,b1,c1,a,bsq,logg,verbos& enddo if (vb) print *,'hhn used mem',2*(nsafe+nintervals) deallocate(mw,mp,poly,poly1) !deallocation may be necessary because these are pointers. + contains + elemental real function w(x) + implicit none + real,intent(in) :: x + w=x**alpha*exp(-sign*x*x) ! inline function for weight + end function w + elemental real function sechb(x) + implicit none + real,intent(in) :: x + sechb=2.*x/(1.+x*x) + end function sechb + elemental real function sech(x) + implicit none + real,intent(in) :: x + sech=sechb(exp(-x)) + end function sech + elemental real function sech2(x) + implicit none + real,intent(in) :: x + sech2=sech(abs(x))**2 + end function sech2 + elemental real function tanhb(x) + implicit none + real,intent(in) :: x + tanhb=2*x/(1+x) + end function tanhb + elemental real function tanhp1(x) + implicit none + real,intent(in) :: x + tanhp1=tanhb(exp(2*x)) ! tanh(x)+1 for x<0 approx. 2*exp(2*x) for x<<-1 + end function tanhp1 + elemental real function tanhm1(x) + implicit none + real,intent(in) :: x + tanhm1=-tanhb(exp(-2*x)) ! tanh(x)-1 for x>0 + end function tanhm1 end subroutine half_hermite_norm subroutine gauss_nodes_functions(n,a,bsq,sp,ortho_f)