istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

added branch for multiple passes (from nocommon rev142); fixed initial condition for second pass, new dispersion module for low density

Browse Source
This commit is contained in:
Daniele Micheletti 2016-01-26 09:14:47 +00:00
parent e1ba175efb
commit c0c966af96
4 changed files with 104 additions and 85 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/beamdata.f90
View File

@ -2,7 +2,7 @@ module beamdata
use const_and_precisions, only : wp_ use const_and_precisions, only : wp_
implicit none implicit none
integer, parameter :: jmx=31,kmx=36,nmx=8000 integer, parameter :: jmx=31,kmx=36,nmx=100000
integer, save :: nrayr,nrayth,nstep integer, save :: nrayr,nrayth,nstep
real(wp_), dimension(:,:,:), allocatable, save :: psjki,ppabs,ccci,tauv,alphav real(wp_), dimension(:,:,:), allocatable, save :: psjki,ppabs,ccci,tauv,alphav
real(wp_), dimension(:,:,:), allocatable, save :: pdjki,currj,didst real(wp_), dimension(:,:,:), allocatable, save :: pdjki,currj,didst

4
src/const_and_precisions.f90
View File

@ -32,7 +32,9 @@
!======================================================================== !========================================================================
integer, parameter :: izero = 0 integer, parameter :: izero = 0
REAL(wp_), PARAMETER :: zero = 0.0_wp_ REAL(wp_), PARAMETER :: zero = 0.0_wp_
REAL(wp_), PARAMETER :: half = 0.5_wp_
REAL(wp_), PARAMETER :: one = 1.0_wp_ REAL(wp_), PARAMETER :: one = 1.0_wp_
REAL(wp_), PARAMETER :: two = 2.0_wp_
real(wp_), parameter :: pi = 3.141592653589793_wp_ ! 3.141592653589793238462643383280 real(wp_), parameter :: pi = 3.141592653589793_wp_ ! 3.141592653589793238462643383280
real(wp_), parameter :: pihalf = 1.57079632679489661923_wp_ real(wp_), parameter :: pihalf = 1.57079632679489661923_wp_
REAL(wp_), PARAMETER :: sqrt_pi = 1.772453850905516_wp_ REAL(wp_), PARAMETER :: sqrt_pi = 1.772453850905516_wp_
@ -102,7 +104,7 @@
REAL(wp_), PARAMETER :: mu0inv = 1._wp_/mu0_ ! REAL(wp_), PARAMETER :: mu0inv = 1._wp_/mu0_ !
! REAL(wp_), PARAMETER :: mc_ = me_*c_ ! [kg*m/s] ! REAL(wp_), PARAMETER :: mc_ = me_*c_ ! [kg*m/s]
! ! f_ce = fce1_*B (B in Tesla): ! ! ! f_ce = fce1_*B (B in Tesla): !
! REAL(wp_), PARAMETER :: wce1_ = e_/me_ ! [rad/s] REAL(wp_), PARAMETER :: wce1_ = e_/me_ ! [rad/s]
! REAL(wp_), PARAMETER :: fce1_ = wce1_/(2*pi) ! [1/s] ! REAL(wp_), PARAMETER :: fce1_ = wce1_/(2*pi) ! [1/s]
! ! f_pl = fpe1_*sqrt(Ne) (Ne in 1/m**3): ! ! ! f_pl = fpe1_*sqrt(Ne) (Ne in 1/m**3): !
! REAL(wp_), PARAMETER :: wpe1_ = 56.4049201 ! [rad/s] ! REAL(wp_), PARAMETER :: wpe1_ = 56.4049201 ! [rad/s]

171
src/dispersion.f90
View File

@ -1,6 +1,7 @@
module dispersion module dispersion
! !
use const_and_precisions, only : wp_,zero,one,im,czero,cunit,pi,sqrt_pi use const_and_precisions, only : wp_,zero,one,half,two,im,czero,cunit,pi, &
sqrt_pi
implicit none implicit none
! local constants ! local constants
integer, parameter :: npts=500 integer, parameter :: npts=500
@ -110,21 +111,28 @@ subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez)
complex(wp_) :: ex,ey,ez,den complex(wp_) :: ex,ey,ez,den
! local variables ! local variables
integer :: i,j,k,imxx,ilrm integer :: i,j,k,imxx,ilrm
real(wp_) :: errnpr,npl2,s,ez2,ey2,ex2,enx2 real(wp_) :: errnpr,npl2,s,ez2,ey2,ex2,enx2,dnl
complex(wp_) :: cc0,cc2,cc4,discr,npra2,npra,npr,npr2,e330, & complex(wp_) :: cc0,cc2,cc4,discr,rdiscr,npra2,npra,npr,npr2,e330,e11,e22, &
e11,e22,e12,e13,e23,a13,a31,a23,a32,a33 e12,e13,e23,a11,a22,a33,a12,a13,a23,a330,a1122,a123,a330n, &
cc0t,cc2t,cc4t
complex(wp_) :: epsl(3,3,lrm),sepsl(3,3) complex(wp_) :: epsl(3,3,lrm),sepsl(3,3)
! !
err=0 err=0
errnpr=one errnpr=one
npra2=nprf**2 npra2=nprf**2
npl2=npl*npl npl2=npl*npl
dnl=one-npl2
!
if(sox>0.and.xg>0.15.and.xg<0.4) then
imxx=1
else
imxx=abs(imx) imxx=abs(imx)
end if
! !
if (fast.eq.1) then if (fast.eq.1) then
call diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm) call diel_tens_wr(yg,mu,npl,a330,epsl,lrm)
else else
call diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast) call diel_tens_fr(yg,mu,npl,a330,epsl,lrm,fast)
end if end if
! !
do do
@ -141,40 +149,64 @@ subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez)
! !
npra=sqrt(npra2) npra=sqrt(npra2)
! !
e11=sepsl(1,1) a11 = xg*sepsl(1,1)
e22=sepsl(2,2) a22 = xg*sepsl(2,2)
e12=sepsl(1,2) a12 = xg*sepsl(1,2)
a33=sepsl(3,3) a33 = xg*sepsl(3,3)
a13=sepsl(1,3) a13 = xg*sepsl(1,3)
a23=sepsl(2,3) a23 = xg*sepsl(2,3)
a31=a13 a330= xg*a330
a32=-a23 ! a31 = a13
! e33=e330+npra2*a33 ! a32 =-a23
e13=npra*a13 !
e23=npra*a23 e11 = one + a11
! e21=-e12 e22 = one + a22
! e31=e13 e12 = a12
! e32=-e23 e330= one + a330
! e33 = e330 + npra2*a33
e13 = npra*a13
e23 = npra*a23
! e21 =-e12
! e31 = e13
! e32 =-e23
! !
! if(i.gt.2.and.errnpr.lt.1.0e-3_wp_) exit ! if(i.gt.2.and.errnpr.lt.1.0e-3_wp_) exit
! !
cc4=(e11-npl2)*(one-a33)+(a13+npl)*(a31+npl) cc4 = (a11 + dnl)*(one - a33) + (a13 + npl)*(a13 + npl)
cc2=-e12*e12*(one-a33)-a32*e12*(a13+npl)+a23*e12*(a31+npl) & cc2 =-a12*a12*(one - a33) + two*a23*a12*(a13 + npl) &
-(a23*a32+e330+(e22-npl2)*(one-a33))*(e11-npl2) & -(e330 + (a22 + dnl)*(one - a33) - a23*a23)*(a11 + dnl) &
-(a13+npl)*(a31+npl)*(e22-npl2) -(a13 + npl)*(a13 + npl)*(a22 + dnl)
cc0=e330*((e11-npl2)*(e22-npl2)+e12*e12) cc0 = e330*((a11 + dnl)*(a22 + dnl) + a12*a12)
! !
discr=cc2*cc2-4.0_wp_*cc0*cc4 cc4t = cc4 - one
cc2t = half*cc2 + dnl
cc0t = cc0 - dnl*dnl
! !
if(yg.gt.one) then a1122 = a11*a22 + a12*a12
a330n = a330 + dnl*a33
a123 = a12*a23 - a13*a22
!
discr = (cc2t*cc2t - cc0t*cc4t) &
-( npl2*a1122 + dnl*a22*a330n + (dnl*a23)**2 + two*dnl*npl*a123 &
+ a1122*a330n + dnl*a13*a123 + dnl*a23*(a12*a13 + a11*a23) )
!
! if(yg.gt.one) then
! s=sox
! if(dimag(discr).LE.zero) s=-s
! else
! s=-sox
! if(dimag(discr).ge.zero) s=-s
! end if
!
rdiscr=sqrt(discr)
if(dimag(rdiscr/cc4).gt.0.0d0) then
! if(dimag(discr).gt.0.0d0) then
s=sox s=sox
if(dimag(discr).LE.zero) s=-s
else else
s=-sox s=-sox
if(dble(discr).le.zero.and.dimag(discr).ge.zero) s=-s
end if end if
! !
npr2=(-cc2+s*sqrt(discr))/(2.0_wp_*cc4) npr2=(s*rdiscr - half*cc2)/cc4
! !
errnpr=abs(one-abs(npr2)/abs(npra2)) errnpr=abs(one-abs(npr2)/abs(npra2))
if(i.gt.1.and.errnpr.lt.1.0e-5_wp_) exit if(i.gt.1.and.errnpr.lt.1.0e-5_wp_) exit
@ -183,35 +215,23 @@ subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez)
if(i.gt.imxx.and.imxx.gt.1) then if(i.gt.imxx.and.imxx.gt.1) then
if (imx.lt.0) then if (imx.lt.0) then
write(*,"(' X =',f7.4,' Y =',f10.7,' Nperp =',f7.4,': convergence & err=1
&disabled.',e12.5)") xg,yg,sqrt(abs(npr2)),npl
imxx=1 imxx=1
else else
write(*,"(' X =',f7.4,' Y =',f10.7,' Nperp =',f7.4,': convergence & err=2
&failed.',e12.5)") xg,yg,sqrt(abs(npr2)),npl
exit exit
end if end if
else else
exit exit
end if end if
print*,yg,imx,imxx
end do end do
! !
! if(i.gt.imx) print*,' i>imx ',yg,errnpr,i if(dble(sqrt(npr2)).lt.zero.or.npr2.ne.npr2.or.abs(npr2).ge.huge(one).or. &
!
if(sqrt(dble(npr2)).lt.zero.or.npr2.ne.npr2.or.abs(npr2).eq.huge(one).or. &
abs(npr2).le.tiny(one)) then abs(npr2).le.tiny(one)) then
write(*,"(' X =',f7.4,' Y =',f7.4,' Nperp =',f7.4,'!')") xg,yg,sqrt(abs(npr2))
npr2=czero npr2=czero
err=err+4
end if end if
! if(dble(npr2).lt.zero) then
! npr2=zero
! print*,' Y =',yg,' npr2 < 0'
! err=99
! end if
!
! write(11,99) yg,dble(npr2),dimag(npr2),nprf**2,dble(i)
! !
npr=sqrt(npr2) npr=sqrt(npr2)
nprr=dble(npr) nprr=dble(npr)
@ -250,7 +270,7 @@ end subroutine warmdisp
! !
! !
! !
subroutine diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast) subroutine diel_tens_fr(yg,mu,npl,a330,epsl,lrm,fast)
! Fully relativistic case computation of dielectric tensor elements ! Fully relativistic case computation of dielectric tensor elements
! up to third order in Larmor radius for hermitian part ! up to third order in Larmor radius for hermitian part
! !
@ -258,8 +278,8 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast)
implicit none implicit none
! arguments ! arguments
integer :: lrm,fast integer :: lrm,fast
real(wp_) :: xg,yg,mu,npl real(wp_) :: yg,mu,npl
complex(wp_) :: e330 complex(wp_) :: a330
complex(wp_), dimension(3,3,lrm) :: epsl complex(wp_), dimension(3,3,lrm) :: epsl
! local variables ! local variables
integer :: i,j,l,is,k,lm integer :: i,j,l,is,k,lm
@ -284,18 +304,11 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast)
end do end do
end do end do
! !
select case(fast) if (fast<4) then
case(2:3)
call hermitian(rr,yg,mu,npl,cr,fast,lrm) call hermitian(rr,yg,mu,npl,cr,fast,lrm)
else
case(4:)
call hermitian_2(rr,yg,mu,npl,cr,fast,lrm) call hermitian_2(rr,yg,mu,npl,cr,fast,lrm)
end if
case default
write(*,*) "unexpected value for flag 'fast' in dispersion:", fast
end select
! !
call antihermitian(ri,yg,mu,npl,ci,lrm) call antihermitian(ri,yg,mu,npl,ci,lrm)
! !
@ -336,19 +349,16 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast)
ca23=ca23+l*asl*cq1p/yg ca23=ca23+l*asl*cq1p/yg
ca33=ca33+asl*cq2p/yg**2 ca33=ca33+asl*cq2p/yg**2
end do end do
epsl(1,1,l) = - xg*ca11*fal epsl(1,1,l) = - ca11*fal
epsl(1,2,l) = + im*xg*ca12*fal epsl(1,2,l) = + im*ca12*fal
epsl(2,2,l) = - xg*ca22*fal epsl(2,2,l) = - ca22*fal
epsl(1,3,l) = - xg*ca13*fal epsl(1,3,l) = - ca13*fal
epsl(2,3,l) = - im*xg*ca23*fal epsl(2,3,l) = - im*ca23*fal
epsl(3,3,l) = - xg*ca33*fal epsl(3,3,l) = - ca33*fal
end do end do
! !
cq2p=rr(0,2,0) cq2p=rr(0,2,0)
e330=one+xg*cq2p a330=cq2p
!
epsl(1,1,1) = one + epsl(1,1,1)
epsl(2,2,1) = one + epsl(2,2,1)
! !
do l=1,lrm do l=1,lrm
epsl(2,1,l) = - epsl(1,2,l) epsl(2,1,l) = - epsl(1,2,l)
@ -915,7 +925,7 @@ end subroutine expinit
! !
! !
! !
subroutine diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm) subroutine diel_tens_wr(yg,mu,npl,a330,epsl,lrm)
! Weakly relativistic dielectric tensor computation of dielectric ! Weakly relativistic dielectric tensor computation of dielectric
! tensor elements (Krivenki and Orefice, JPP 30,125 - 1983) ! tensor elements (Krivenki and Orefice, JPP 30,125 - 1983)
! !
@ -923,8 +933,8 @@ subroutine diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm)
implicit none implicit none
! arguments ! arguments
integer :: lrm integer :: lrm
real(wp_) :: xg,yg,npl,mu real(wp_) :: yg,npl,mu
complex(wp_) :: e330,epsl(3,3,lrm) complex(wp_) :: a330,epsl(3,3,lrm)
! local variables ! local variables
integer :: l,lm,is,k integer :: l,lm,is,k
real(wp_) :: npl2,fcl,w,asl,bsl real(wp_) :: npl2,fcl,w,asl,bsl
@ -964,19 +974,16 @@ subroutine diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm)
ca23=ca23+l*asl*cq1p/yg ca23=ca23+l*asl*cq1p/yg
ca33=ca33+asl*cq2p/yg**2 ca33=ca33+asl*cq2p/yg**2
end do end do
epsl(1,1,l) = -xg*ca11*fcl epsl(1,1,l) = -ca11*fcl
epsl(1,2,l) = +im*xg*ca12*fcl epsl(1,2,l) = +im*ca12*fcl
epsl(2,2,l) = -xg*ca22*fcl epsl(2,2,l) = -ca22*fcl
epsl(1,3,l) = -xg*ca13*fcl epsl(1,3,l) = -ca13*fcl
epsl(2,3,l) = -im*xg*ca23*fcl epsl(2,3,l) = -im*ca23*fcl
epsl(3,3,l) = -xg*ca33*fcl epsl(3,3,l) = -ca33*fcl
end do end do
! !
cq2p=cefp(0,1)+mu*npl2*(cefp(0,2)+cefp(0,0)-2.0_wp_*cefp(0,1)) cq2p=cefp(0,1)+mu*npl2*(cefp(0,2)+cefp(0,0)-2.0_wp_*cefp(0,1))
e330=one-xg*mu*cq2p a330=-mu*cq2p
!
epsl(1,1,1) = one + epsl(1,1,1)
epsl(2,2,1) = one + epsl(2,2,1)
! !
do l=1,lrm do l=1,lrm
epsl(2,1,l) = - epsl(1,2,l) epsl(2,1,l) = - epsl(1,2,l)

12
src/gray.f
View File

@ -163,6 +163,7 @@ c
write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabstot write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabstot
currtka =currtot*1.0e3_wp_ currtka =currtot*1.0e3_wp_
write(*,'(a,f9.4)') 'I_tot (kA) = ',currtka write(*,'(a,f9.4)') 'I_tot (kA) = ',currtka
write(7,99) currtka,pabstot
c c
if (index_rt.eq.1) then if (index_rt.eq.1) then
if(iequil.eq.2) write(*,*) 'EQUILIBRIUM CASE : ',filenmeqq if(iequil.eq.2) write(*,*) 'EQUILIBRIUM CASE : ',filenmeqq
@ -327,12 +328,16 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
end if end if
call wall_refl(xv-dst*anv,anv,ext(j,k,iop(j,k)), call wall_refl(xv-dst*anv,anv,ext(j,k,iop(j,k)),
. eyt(j,k,iop(j,k)),xvrfl,anvrfl,extr,eytr,anw,irfl) . eyt(j,k,iop(j,k)),xvrfl,anvrfl,extr,eytr,anw,irfl)
if(index_rt.eq.1) then
istore(j,k)=istore(j,k)+1 istore(j,k)=istore(j,k)+1
yyrfl(j,k,1:3)=xvrfl yyrfl(j,k,1:3)=xvrfl
yyrfl(j,k,4:6)=anvrfl yyrfl(j,k,4:6)=anvrfl
tau1v(j,k)=tauv(j,k,iiv(j,k)) tau1v(j,k)=tauv(j,k,iiv(j,k))
ext(j,k,iop(j,k))=extr ext(j,k,iop(j,k))=extr
eyt(j,k,iop(j,k))=eytr eyt(j,k,iop(j,k))=eytr
end if
if (j.lt.jclosest) then if (j.lt.jclosest) then
jclosest=j jclosest=j
anwcl=anw anwcl=anw
@ -3113,6 +3118,7 @@ c
common/anv/anv common/anv/anv
common/xv/xv common/xv/xv
c c
ierr=0
xx=y(1) xx=y(1)
yy=y(2) yy=y(2)
zz=y(3) zz=y(3)
@ -4410,7 +4416,7 @@ c
alpha=2.0_wp_*akim*ratiovgr alpha=2.0_wp_*akim*ratiovgr
if(alpha.lt.0.0_wp_) then if(alpha.lt.0.0_wp_) then
ierr=94 ierr=94
print*,' IERR = ', ierr,' alpha negative' print*,' IERR = ', ierr,' alpha negative', alpha
end if end if
c c
if(ieccd.gt.0) then if(ieccd.gt.0) then
@ -5175,6 +5181,7 @@ c
use const_and_precisions, only : wp_ use const_and_precisions, only : wp_
use reflections, only : inters_linewall,inside,rlim,zlim,nlim use reflections, only : inters_linewall,inside,rlim,zlim,nlim
use graydata_flags, only : dst use graydata_flags, only : dst
use beamdata, only : nstep
implicit none implicit none
c arguments c arguments
integer :: ivac integer :: ivac
@ -5222,6 +5229,9 @@ c used by fwork!!!
smax=smax*1.0e2_wp_+st smax=smax*1.0e2_wp_+st
end if end if
i=0 i=0
if(smax.gt.dst*nstep) then
smax=dst*nstep
end if
do do
st=i*dst st=i*dst
xvend=xv0+st*anv0 xvend=xv0+st*anv0