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

trunk/graycore: minor fixes

Browse Source
This commit is contained in:
Daniele Micheletti 2019-03-28 09:50:28 +00:00
parent 27f1793f14
commit f600108884
2 changed files with 67 additions and 62 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

63
src/graycore.f90
View File

@ -185,7 +185,7 @@ contains
nbeam_pass=1 ! max n of beam per pass
index_rt=0 ! global beam index: 1,O 2,X 1st pass
! / \ / \
! | | | |
do ip=1,ipass ! 3,O 4,X 5,O 6,X 2nd pass
pabs_pass = zero
@ -194,7 +194,7 @@ contains
nbeam_pass = 2*nbeam_pass ! max n of beams in current pass
if(ip.gt.1) then
du1 = zero ! igrad=0 from 2nd pass
du1 = zero
gri = zero
ggri = zero
if(ip.eq.ipass) cpl = (/zero,zero/) ! no successive passes
@ -220,12 +220,12 @@ contains
call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
call print_headers((/' '/),index_rt)
if(ip.eq.1) then ! first pass
igrad_b = igrad
if(ip.eq.1) then ! 1st pass
igrad_b = igrad ! * input value, igrad_b=0 from 2nd pass
tau1 = zero
tau1 = zero ! * tau from previous passes
etau1 = one
cpl1 = one
cpl1 = one ! * coupling from previous passes
lgcpl1 = zero
p0ray = p0jk ! * initial beam power
@ -233,7 +233,6 @@ contains
call set_pol(yw,bres,sox,psipol,chipol,ext,eyt) ! * initial polarization
do jk=1,nray
if(iwait(jk)) cycle
zzm = yw(3,jk)*0.01_wp_
rrm = sqrt(yw(1,jk)*yw(1,jk)+yw(2,jk)*yw(2,jk))*0.01_wp_
@ -244,7 +243,7 @@ contains
end if
end do
else ! successive passes
else ! 2nd+ passes
ipar = (index_rt+1)/2-1 ! * parent beam index
yw = yynext(:,:,ipar) ! * starting coordinates from
ypw = yypnext(:,:,ipar) ! parent beam last step
@ -307,7 +306,7 @@ contains
if(ent_pl) then ! ray enters plasma
call plasma_in(jk,xv,anv,bres,sox,cpl,psipol,chipol,iop,ext,eyt)
if(iop(jk).eq.1 .and. ip==1) then ! * 1st entrance on 1st pass (ray never entered in plasma) => continue current pass
if(iop(jk).eq.1 .and. ip==1) then ! * 1st entrance on 1st pass (ray hasn't entered in plasma yet) => continue current pass
if(ipol.eq.0) then ! + IF single mode propagation
cpl = cpl0
@ -333,29 +332,29 @@ contains
igrad_b = 0 ! + switch to ray-tracing
iwait(jk) = .true. ! + stop advancement and H&CD computation for current ray
if(ip.lt.ipass) then ! not last pass
yynext(:,jk,index_rt) = yw0(:,jk) ! + copy starting coordinates
if(ip.lt.ipass) then ! + not last pass
yynext(:,jk,index_rt) = yw0(:,jk) ! . copy starting coordinates
yypnext(:,jk,index_rt) = ypw0(:,jk) ! for next pass from last step
stnext(jk,index_rt) = stv(jk) - dst ! + starting step for next pass = last step
stnext(jk,index_rt) = stv(jk) - dst ! . starting step for next pass = last step
if(cpl(1).lt.etaucr) then ! + low coupled power for O-mode => de-activate derived rays
if(cpl(1).lt.etaucr) then ! . low coupled power for O-mode => de-activate derived rays
call turnoffray(jk,ip+1,2*ib-1,iroff)
if(cpl(1).le.comp_tiny) cpl(1)=etaucr
end if
if(cpl(2).lt.etaucr) then ! + low coupled power for X-mode => de-activate derived rays
if(cpl(2).lt.etaucr) then ! . low coupled power for X-mode => de-activate derived rays
call turnoffray(jk,ip+1,2*ib,iroff)
if(cpl(2).le.comp_tiny) cpl(2)=etaucr
end if
taus(jk,iO:iO+1) = tau1(jk) + tau0(jk) ! + starting tau for next O-mode pass
cpls(jk,iO) = cpl1(jk) * cpl(1) ! + cumulative coupling for next O-mode pass
cpls(jk,iO+1) = cpl1(jk) * cpl(2) ! + cumulative coupling for next X-mode pass
taus(jk,iO:iO+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass
cpls(jk,iO) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-mode pass
cpls(jk,iO+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass
if(jk.eq.1) then ! + polarization angles at plasma boundary for central ray
if(jk.eq.1) then ! . polarization angles at plasma boundary for central ray
psipv(iO:iO+1) = psipol
chipv(iO:iO+1) = chipol
end if
else
else ! * 1st entrance on 2nd+ pass (ray hasn't entered in plasma since end of previous pass) => continue current pass
cpl = (/zero,zero/)
end if
end if
@ -365,7 +364,7 @@ contains
end if
if(ent_wl) then ! ray enters vessel
iow(jk)=iow(jk)+1
iow(jk)=iow(jk)+1 ! * out->in
else if(ext_wl) then ! ray exit vessel
call wall_out(jk,ins_pl,xv,anv,bres,sox,psipol,chipol,iow,iop,ext,eyt)
@ -385,7 +384,7 @@ contains
chipv(index_rt) = chipol
end if
if(ins_pl) then ! * plasma-wall overlapping => wall crossing=plasma crossing => end of current pass
if(ins_pl) then ! * plasma-wall overlapping => wall+plasma crossing => end of current pass
iwait(jk) = .true. ! + stop advancement and H&CD computation for current ray
if(ip.lt.ipass) then ! + not last pass
@ -418,11 +417,11 @@ contains
iopmin = min(iopmin,iop(jk))
if(ip.lt.ipass) then ! not last pass
yw0(:,jk) = yw(:,jk) ! * store current coordinates
ypw0(:,jk) = ypw(:,jk) ! in case pass ends on next step
yw0(:,jk) = yw(:,jk) ! * store current coordinates in case
ypw0(:,jk) = ypw(:,jk) ! current pass ends on next step
end if
! compute ECRH&CD
! compute ECRH&CD if (inside plasma & power available>0 & ray still active)
if(ierrn==0 .and. iwarm>0 .and. ins_pl .and. &
(tau1(jk)+tau0(jk)+lgcpl1(jk))<=taucr .and. .not.iwait(jk)) then ! H&CD computation check
tekev=temp(psinv)
@ -456,14 +455,14 @@ contains
dids0(jk)=dids
ccci0(jk)=ccci(jk,i)
if(iwait(jk)) then
if(iwait(jk)) then ! copy values from last pass for inactive ray
ppabs(jk,i:nstep) = ppabs(jk,i-1)
ccci(jk,i:nstep) = ccci(jk,i-1)
psjki(jk,i:nstep) = psjki(jk,i-1)
else
call print_output(i,jk,stv(jk),p0ray(jk)/etau1(jk),xv,psinv, &
btot,bv,ak0,anpl,anpr,anv,anprim,dens,tekev,alpha,tau,dids, &
nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg) !****** p0ray(jk)/etau1(jk) fattore normalizzazione per dids da controllare
nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg) ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes)
end if
end do
@ -490,7 +489,7 @@ contains
istop_pass = istop_pass +1 ! * +1 non propagating beam
if(ip.lt.ipass) call turnoffray(0,ip,ib,iroff) ! * de-activate derived beams
exit
else if(all(iwait)) then ! all rays in current beam are waiting for next pass
else if(all(iwait)) then ! all rays in current beam are waiting for next pass => do not increase istop_pass
exit
end if
end do
@ -514,10 +513,10 @@ contains
icd_pass(iox) = icd_pass(iox) + icd_beam
if(ip.lt.ipass) then
cpl_beam1 = sum(p0ray * exp(-tau0) * cpls(:,iO)/cpl1) / sum(p0ray * exp(-tau0)) ! average O-mode coupling
cpl_beam2 = one - cpl_beam1
cpl_cbeam1 = cpls(1,iO)/cpl1(1)
cpl_cbeam2 = one - cpl_cbeam1
cpl_beam1 = sum(p0ray * exp(-tau0) * cpls(:,iO)/cpl1) / sum(p0ray * exp(-tau0)) ! average O-mode coupling for next beam
cpl_beam2 = one - cpl_beam1 ! average X-mode coupling for next beam
cpl_cbeam1 = cpls(1,iO)/cpl1(1) ! central ray O-mode coupling for next beam
cpl_cbeam2 = one - cpl_cbeam1 ! central ray X-mode coupling for next beam
else
cpl_beam1 = zero
cpl_beam2 = zero
@ -533,7 +532,7 @@ contains
write(*,'(a,f9.4)') 'I_tot (kA) = ',icd_beam*1.0e3_wp_
if(ip.lt.ipass) then
write(*,'(a,2(f9.4,1x))') 'Coupling (average, O/X):',cpl_beam1,cpl_beam2 ! average coupling for next O/X beams
write(*,'(a,2(f9.4,1x))') 'Coupling (c. beam, O/X):',cpl_cbeam1,cpl_cbeam2 ! central beam coupling for next O/X beams
write(*,'(a,2(f9.4,1x))') 'Coupling (ctr ray, O/X):',cpl_cbeam1,cpl_cbeam2 ! central ray coupling for next O/X beams
end if
call print_pec(rhop_tab,rhot_tab,jphi_beam,jcd_beam,dpdv_beam,currins_beam, &

64
src/multipass.f90
View File

@ -8,28 +8,29 @@ module multipass
implicit none
integer, save :: nbeam_max
integer, save :: nbeam_tot
integer, save :: nbeam_max ! max n of beams active at a time
integer, save :: nbeam_tot ! total n of beams
contains
! ------------------------------
subroutine plasma_in(i,xv,anv,bres,sox,cpl,psipol1,chipol1,iop,ext,eyt) ! ray enters plasma
implicit none
! arguments
integer, intent(in) :: i ! ray index
real(wp_), dimension(3), intent(in) :: xv,anv
real(wp_), intent(in) :: bres,sox
real(wp_), dimension(2), intent(out) :: cpl ! coupling (O/X)
real(wp_), intent(out) :: psipol1,chipol1
integer, dimension(:), intent(inout), pointer :: iop
integer, dimension(:), intent(inout), pointer :: iop ! in/out plasma flag
complex(wp_), dimension(:), intent(inout), pointer :: ext,eyt
!
! local variables
real(wp_) :: rm,csphi,snphi,bphi,br,bz
real(wp_) :: qq1,uu1,vv1,qq,uu,vv,powcpl
real(wp_), dimension(3) :: bv,xmv
complex(wp_) :: ext1,eyt1
!
iop(i)=iop(i)+1
iop(i)=iop(i)+1 ! out->in
xmv=xv*0.01_wp_ ! convert from cm to m
rm=sqrt(xmv(1)**2+xmv(2)**2)
@ -43,7 +44,7 @@ contains
call pol_limit(anv,bv,bres,sox,ext1,eyt1)
call stokes_ce(ext1,eyt1,qq1,uu1,vv1) ! stokes parameter at plasma entrance
call stokes_ce(ext(i),eyt(i),qq,uu,vv) ! stokes parameter at plasma exit/wall reflection
powcpl = half*(one + vv*vv1+uu*uu1+qq*qq1)
powcpl = half*(one + vv*vv1+uu*uu1+qq*qq1) ! coupling for incoming mode
if(sox.eq.-one) then ! incoming mode = O
cpl=(/powcpl,one-powcpl/)
@ -51,7 +52,7 @@ contains
cpl=(/one-powcpl,powcpl/)
end if
if(i.eq.1) then
if(i.eq.1) then ! polarization angles at plasma entrace for central ray
call polellipse(qq1,uu1,vv1,psipol1,chipol1)
psipol1=psipol1/degree ! convert from rad to degree
chipol1=chipol1/degree
@ -63,16 +64,17 @@ contains
! ------------------------------
subroutine plasma_out(i,xv,anv,bres,sox,iop,ext,eyt) ! ray exits plasma
implicit none
! arguments
integer, intent(in) :: i ! ray index
real(wp_), dimension(3), intent(in) :: xv,anv
real(wp_), intent(in) :: bres,sox
integer, dimension(:), intent(inout), pointer :: iop
integer, dimension(:), intent(inout), pointer :: iop ! in/out plasma flag
complex(wp_), dimension(:), intent(out), pointer :: ext,eyt
!
! local variables
real(wp_) :: rm,csphi,snphi,bphi,br,bz
real(wp_), dimension(3) :: bv,xmv
!
iop(i)=iop(i)+1
iop(i)=iop(i)+1 ! in->out
xmv=xv*0.01_wp_ ! convert from cm to m
rm=sqrt(xmv(1)**2+xmv(2)**2)
@ -95,23 +97,24 @@ contains
! ------------------------------
subroutine wall_out(i,ins,xv,anv,bres,sox,psipol1,chipol1,iow,iop,ext,eyt) ! ray exits vessel
implicit none
! arguments
integer, intent(in) :: i ! ray index
logical, intent(in) :: ins ! inside plasma? (ins=1 plasma/wall overlap)
real(wp_), dimension(3), intent(inout) :: xv,anv
real(wp_), intent(in) :: bres,sox
real(wp_), intent(out) :: psipol1,chipol1
integer, dimension(:), intent(inout), pointer :: iow,iop
integer, dimension(:), intent(inout), pointer :: iow,iop ! in/out vessel and plasma flags
complex(wp_), dimension(:), intent(inout), pointer :: ext,eyt
!
! local variables
integer :: irfl
real(wp_) :: qq1,uu1,vv1
real(wp_), dimension(3) :: xvrfl,anvrfl,walln
complex(wp_) :: ext1,eyt1
!
iow(i)=iow(i)+1
if(ins) call plasma_out(i,xv,anv,bres,sox,iop,ext,eyt)
call wall_refl(xv-dst*anv,anv,ext(i),eyt(i),xvrfl,anvrfl,ext1,eyt1,walln,irfl)
ext(i) = ext1
!
iow(i)=iow(i)+1 ! out->in
if(ins) call plasma_out(i,xv,anv,bres,sox,iop,ext,eyt) ! plasma-wall overlapping
call wall_refl(xv-dst*anv,anv,ext(i),eyt(i),xvrfl,anvrfl,ext1,eyt1,walln,irfl) ! ray reflects at wall
ext(i) = ext1 ! save parameters at wall reflection
eyt(i) = eyt1
xv = xvrfl
anv = anvrfl
@ -131,8 +134,9 @@ contains
subroutine initbeam(i,iroff,iboff,iwait,stv,jphi_beam,pins_beam,currins_beam, &
dpdv_beam,jcd_beam) ! initialization at beam propagation start
implicit none
! arguments
integer, intent(in) :: i ! beam index
logical, dimension(:,:), intent(in), pointer :: iroff
logical, dimension(:,:), intent(in), pointer :: iroff ! global ray status (F = active, T = inactive)
logical, intent(out) :: iboff
logical, dimension(:), intent(out), pointer :: iwait
real(wp_), dimension(:), intent(out), pointer :: jphi_beam,pins_beam, &
@ -144,7 +148,7 @@ contains
iboff = .true.
write(*,*)
write(*,'("Beam ",i5," inactive")'),i
else if(.not.all(.not.iwait)) then ! no all rays active
else if(.not.all(.not.iwait)) then ! only some rays active
write(*,*)
write(*,'("WARNING: not all rays in beam ",i5," are active")'),i
end if
@ -161,25 +165,26 @@ contains
subroutine initmultipass(i,iox,iroff,yynext,yypnext,yw0,ypw0,stnext,p0ray, &
taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv) ! initialization before pass loop
implicit none
! arguments
integer, intent(in) :: i, iox ! ipol, mode active on 1st pass
logical, dimension(:,:), intent(out), pointer :: iroff
logical, dimension(:,:), intent(out), pointer :: iroff ! global ray status (F = active, T = inactive)
real(wp_), dimension(:), intent(out), pointer :: p0ray,tau1,etau1,cpl1,lgcpl1, &
psipv,chipv
real(wp_), dimension(:,:), intent(out), pointer :: yw0,ypw0,stnext,taus,cpls
real(wp_), dimension(:,:,:), intent(out), pointer :: yynext,yypnext
!
iroff = .false. ! global ray status (F = active, T = inactive)
if(i.eq.0) call turnoffray(0,1,3-iox,iroff) ! ipol = 0 => stop other mode (iox=1/2 -> stop ib=2/1)
if(i.eq.0) call turnoffray(0,1,3-iox,iroff) ! ipol = 0 => stop other mode (iox=1/2 -> stop ib=2/1 at first pass)
yynext = zero ! starting beam coordinates (1)
yypnext = zero ! starting beam coordinates (2)
yw0 = zero ! temporary beam coordinates (1)
ypw0 = zero ! temporary beam coordinates (2)
stnext = zero ! starting beam step
p0ray = zero ! starting beam power
taus = zero ! starting beam tau
tau1 = zero ! total beam tau
taus = zero ! beam tau from previous passes
tau1 = zero
etau1 = one
cpls = one
cpls = one ! beam coupling from previous passes
cpl1 = one
lgcpl1 = zero
psipv = zero ! psi polarization angle at vacuum-plasma boundary
@ -188,15 +193,16 @@ contains
! ------------------------------
subroutine turnoffray(jk,ip,ib,iroff) ! turn off ray propagation
implicit none
! arguments
integer, intent(in) :: jk, ip, ib ! ray (0=all rays), pass, beam indexes
logical, dimension(:,:), intent(out), pointer :: iroff
!
logical, dimension(:,:), intent(out), pointer :: iroff ! global ray status (F = active, T = inactive)
! local variables
integer :: ipx, i1, i2
!
if(jk==0) then ! stop all rays
do ipx=ip,ipass
i1 = 2**ipx-2+2**(ipx-ip)*(ib-1)+1
i2 = 2**ipx-2+2**(ipx-ip)*ib
do ipx=ip,ipass ! from pass ip to last pass
i1 = 2**ipx-2+2**(ipx-ip)*(ib-1)+1 ! first derived beam at pass ipx
i2 = 2**ipx-2+2**(ipx-ip)*ib ! last derived beam at pass ipx (i1=i2 for ipx=ip)
iroff(:,i1:i2) = .true.
end do
else ! only stop ray jk