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

pec subroutine: new algorithm to compute Jphi

Browse Source
This commit is contained in:
Lorenzo Figini 2014-01-27 15:38:17 +00:00
parent d21f9b12f4
commit 1a60867b5b
2 changed files with 141 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
scripts/beta-alpha_beam.sh
View File

@ -1,9 +1,11 @@
#!/bin/bash #!/bin/bash
# set input/working/output folders # set input/working/output folders
in=/Users/daniela/Desktop/scenario_001/files5s in=/Users/lorenzo/ITER/UL-Grant-161/performance-analysis/case006/in
work=/Users/daniela/Desktop/scenario_001/tmpab_beam work=/Users/lorenzo/Desktop/tmp/sc6
out=/Users/daniela/Desktop/scenario_001/outab_beam_550 out=/Users/lorenzo/Public/case006-mappe-ab-beam
scen='006'
# launching mirrors # launching mirrors
mirror=( "USM" "LSM" ) mirror=( "USM" "LSM" )
@ -58,8 +60,8 @@ cp $in/${base}.prf ${base}.prf
#loop on USM/LSM mirrors #loop on USM/LSM mirrors
for i in 0 1; do for i in 0 1; do
# set output file names # set output file names
id="${mirror[i]}_t$t" id="${mirror[i]}_Scen${scen}_t$t"
f7file=$out/f7_${id}.txt f7file=$out/ECCD_vs_ba_${id}.txt
f48file=$out/f48_${id}.txt f48file=$out/f48_${id}.txt
rm -f $out/log_${id}.txt rm -f $out/log_${id}.txt
# add header in output file # add header in output file

152
src/gray.f
View File

@ -512,9 +512,9 @@ c
write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1' write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1'
write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt' write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt'
write(12,*) ' #i sst psi w1 w2' write(12,*) ' #i sst psi w1 w2'
write(7,*)'#Icd Pa Jphimx dPdVmx '// write(7,*)'#Icd Pa Jphip dPdVp '//
.'rhotj rhotjava rhotp rhotpav '// .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav '//
.'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt' .'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx'
write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt' write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt'
else else
@ -1118,6 +1118,7 @@ c
common/cofffp/cfp common/cofffp/cfp
common/fpas/fpolas common/fpas/fpolas
common/rhotsx/rhotsx common/rhotsx/rhotsx
common/phitedge/phitedge
common/rrtor/rrtor common/rrtor/rrtor
common/cnt/rup,zup,rlw,zlw common/cnt/rup,zup,rlw,zlw
common/limiter/rlim,zlim,nlim common/limiter/rlim,zlim,nlim
@ -1457,8 +1458,10 @@ c compute B_toroidal on axis
c compute normalized rho_tor from eqdsk q profile c compute normalized rho_tor from eqdsk q profile
call rhotor(nr) call rhotor(nr)
c phitedge=abs(psia)*rhotsx*2*pi phitedge=abs(psia)*rhotsx*2*pi
c rrtor=sqrt(phitedge/abs(btrcen)/pi) rrtor=sqrt(phitedge/abs(btrcen)/pi)
call rhopol
c print*,rhotsx,phitedge,rrtor,abs(psia)
c compute flux surface averaged quantities c compute flux surface averaged quantities
@ -1822,6 +1825,7 @@ c
common/qpsi/qpsi common/qpsi/qpsi
common/rhotsx/rhotsx common/rhotsx/rhotsx
common/crhot/crhot common/crhot/crhot
common/cq/cq
c c
c normalized toroidal rho : ~ Integral q(psi) dpsi c normalized toroidal rho : ~ Integral q(psi) dpsi
c c
@ -1847,6 +1851,21 @@ c
return return
end end
function fq_eq(psi)
implicit real*8(a-h,o-z)
parameter(nnw=501)
dimension psinr(nnw),cq(nnw,4)
common/psinr/psinr
common/eqnn/nr,nz,npp,nintp
common/cq/cq
irt=int((nr-1)*psi+1)
if(irt.eq.0) irt=1
if(irt.eq.nr) irt=nr-1
dps=psi-psinr(irt)
fq_eq=spli(cq,nr,irt,dps)
return
end
function frhotor_eq(psi) function frhotor_eq(psi)
implicit real*8(a-h,o-z) implicit real*8(a-h,o-z)
parameter(nnw=501) parameter(nnw=501)
@ -1883,7 +1902,57 @@ c
if(ip.eq.nintp) ip=nintp-1 if(ip.eq.nintp) ip=nintp-1
dps=rpsi-rpstab(ip) dps=rpsi-rpstab(ip)
frhotor_av=spli(crhotq,nintp,ip,dps) frhotor_av=spli(crhotq,nintp,ip,dps)
return
end
subroutine rhopol
implicit real*8(a-h,o-z)
parameter(nnr=101,nrest=nnr+4)
parameter(lwrkp=nnr*4+nrest*16)
dimension rhop(nnr),rhot(nnr),rhopi(nnr)
dimension trp(nrest),crp(nrest),wp(nrest)
dimension wrkp(lwrkp),iwrkp(nrest)
common/coffrtp/trp
common/coffrn/nsrp
common/coffrp/crp
dr=1.0d0/dble(nnr-1)
do i=1,nnr
rhop(i)=(i-1)*dr
psin=rhop(i)*rhop(i)
rhot(i)=frhotor(psin)
wp(i)=1.0d0
end do
wp(1)=1.0d3
wp(nnr)=1.0d3
c spline interpolation of rhopol versus rhotor
iopt=0
xb=0.0d0
xe=1.0d0
ss=0.00001
kspl=3
call curfit(iopt,nnr,rhot,rhop,wp,xb,xe,kspl,ss,nrest,nsrp,
. trp,crp,rp,wrkp,lwrkp,iwrkp,ier)
c print*,ier
call splev(trp,nsrp,crp,3,rhot,rhopi,nnr,ier)
do i=1,nnr
write(54,*) rhop(i),rhot(i),rhopi(i)
end do
return
end
function frhopol(rhot)
implicit real*8(a-h,o-z)
parameter(nnr=101,nrest=nnr+4)
dimension trp(nrest),crp(nrest),rrs(1),ffspl(1)
common/coffrtp/trp
common/coffrn/nsrp
common/coffrp/crp
rrs(1)=rhot
call splev(trp,nsrp,crp,3,rrs,ffspl,1,ier)
frhopol=ffspl(1)
return return
end end
@ -2158,6 +2227,8 @@ c
dimension pstab(nnintp),varea(nnintp),vvol(nnintp),rpstab(nnintp) dimension pstab(nnintp),varea(nnintp),vvol(nnintp),rpstab(nnintp)
dimension rri(nnintp),rbav(nnintp),bav(nnintp),rhotqv(nnintp) dimension rri(nnintp),rbav(nnintp),bav(nnintp),rhotqv(nnintp)
dimension dadrhotv(nnintp),cdadrhot(nnintp,4)
dimension dvdrhotv(nnintp),cdvdrhot(nnintp,4)
dimension bmxpsi(nnintp),bmnpsi(nnintp),ffc(nnintp) dimension bmxpsi(nnintp),bmnpsi(nnintp),ffc(nnintp)
dimension vcurrp(nnintp),vajphiav(nnintp),qqv(nnintp) dimension vcurrp(nnintp),vajphiav(nnintp),qqv(nnintp)
dimension rcon(2*ncnt+1),zcon(2*ncnt+1) dimension rcon(2*ncnt+1),zcon(2*ncnt+1)
@ -2186,6 +2257,9 @@ c
common/cnt/rup,zup,rlw,zlw common/cnt/rup,zup,rlw,zlw
common/bound/zbmin,zbmax common/bound/zbmin,zbmax
common/rarea/rarea common/rarea/rarea
common/phitedge/phitedge
common/cdadrhot/cdadrhot
common/cdvdrhot/cdvdrhot
c c
common/coffvl/ch common/coffvl/ch
common/coffdvl/ch01 common/coffdvl/ch01
@ -2368,6 +2442,11 @@ c ratio_pltor = Jcd_||/J_phi Jcd_|| = <Jcd.b>
vratjb(jp)=ratio_cdbtor vratjb(jp)=ratio_cdbtor
qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi)) qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi))
qqv(jp)=qq qqv(jp)=qq
dadrhotv(jp)=phitedge*frhotor_eq(height2)/fq_eq(height2)
. *dadpsi/pi
dvdrhotv(jp)=phitedge*frhotor_eq(height2)/fq_eq(height2)
. *dvdpsi/pi
c c
c write(57,99) sqrt(pstab(jp)),pstab(jp),riav,dvdpsi,area,vvol(jp) c write(57,99) sqrt(pstab(jp)),pstab(jp),riav,dvdpsi,area,vvol(jp)
@ -2422,7 +2501,7 @@ c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/<sqrt(1-lam*B(rhop)/Bmx)>
end do end do
write(56,*)' #psi rhot_eq rhot_av |<B>| |Bmx| |Bmn|'// write(56,*)' #psi rhot_eq rhot_av |<B>| |Bmx| |Bmn|'//
.' Area Vol |I_pl| <J_phi> qq fc ratio_cdbtor' .' Area Vol |I_pl| <J_phi> qq fc'
qqv(1)=qqv(2) qqv(1)=qqv(2)
vajphiav(1)=vajphiav(2) vajphiav(1)=vajphiav(2)
@ -2467,6 +2546,10 @@ c spline coefficients of rhot
call difcs(rpstab,ffc,nintp,iopt,cfc,ier) call difcs(rpstab,ffc,nintp,iopt,cfc,ier)
iopt=0 iopt=0
call difcs(rpstab,rhotqv,nintp,iopt,crhotq,ier) call difcs(rpstab,rhotqv,nintp,iopt,crhotq,ier)
iopt=0
call difcs(rpstab,dadrhotv,nintp,iopt,cdadrhot,ier)
iopt=0
call difcs(rpstab,dvdrhotv,nintp,iopt,cdvdrhot,ier)
c spline interpolation of H(lambda,rhop) and dH/dlambda c spline interpolation of H(lambda,rhop) and dH/dlambda
@ -2485,9 +2568,37 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda
return return
99 format(20(1x,e12.5)) 99 format(20(1x,e12.5))
end end
c
c function fdadrhot(rpsi)
c implicit real*8(a-h,o-z)
parameter(nnintp=101)
dimension rpstab(nnintp),cdadrhot(nnintp,4)
common/pstab/rpstab
common/eqnn/nr,nz,npp,nintp
common/cdadrhot/cdadrhot
ip=int((nintp-1)*rpsi+1)
if(ip.eq.0) ip=1
if(ip.eq.nintp) ip=nintp-1
dps=rpsi-rpstab(ip)
fdadrhot=spli(cdadrhot,nintp,ip,dps)
return
end
function fdvdrhot(rpsi)
implicit real*8(a-h,o-z)
parameter(nnintp=101)
dimension rpstab(nnintp),cdvdrhot(nnintp,4)
common/pstab/rpstab
common/eqnn/nr,nz,npp,nintp
common/cdvdrhot/cdvdrhot
ip=int((nintp-1)*rpsi+1)
if(ip.eq.0) ip=1
if(ip.eq.nintp) ip=nintp-1
dps=rpsi-rpstab(ip)
fdvdrhot=spli(cdvdrhot,nintp,ip,dps)
return
end
subroutine vectinit subroutine vectinit
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(jmx=31,kmx=36,nmx=8000) parameter(jmx=31,kmx=36,nmx=8000)
@ -6180,7 +6291,12 @@ c of gaussian profile
ajcdbv(i)=ajphiv(i)*ratjbv(i) ajcdbv(i)=ajphiv(i)*ratjbv(i)
ajplv(i)=ajphiv(i)*ratjplv(i) ajplv(i)=ajphiv(i)*ratjplv(i)
end do end do
rhpp=frhopol(rhotpav)
rhpj=frhopol(rhotjava)
dpdvp=pabs*2.0d0/(sqrt(pi)*drhotjava*fdvdrhot(rhpp))
ajphip=currt*2.0d0/(sqrt(pi)*drhotjava*fdadrhot(rhpj))
call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx, call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx,
. drhotp,drhopp) . drhotp,drhopp)
if(ieccd.ne.0) then if(ieccd.ne.0) then
@ -6253,18 +6369,18 @@ c dPdV [MW/m^3], Jcd [MA/m^2]
currtka=currt*1.0d3 currtka=currt*1.0d3
write(6,*)' ' write(6,*)' '
write(6,*)'#beta0 alpha0 Icd Pa Jphimx dPdVmx '// write(6,*)'#beta0 alpha0 Icd Pa Jphip dPdVp '//
.'rhotj rhotjava rhotp rhotpav '// .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav '//
.'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt' .'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx'
write(6,99) beta0,alpha0,currtka,pabstot,ajmxfi,dpdvmx, write(6,99) beta0,alpha0,currtka,pabstot,ajphip,dpdvp,
. rhotjfi,rhotjava,rhotp,rhotpav, . rhotjfi,rhotjava,rhotp,rhotpav,
. drhotjava,drhotpav, . drhotjava,drhotpav,
. stmx,psipol,chipol,real(index_rt) . stmx,psipol,chipol,real(index_rt),ajmxfi,dpdvmx
write(7,99) currtka,pabstot,ajmxfi,dpdvmx, write(7,99) currtka,pabstot,ajphip,dpdvp,
. rhotjfi,rhotjava,rhotp,rhotpav, . rhotjfi,rhotjava,rhotp,rhotpav,
. drhotjava,drhotpav,ratjbmx, . drhotjava,drhotpav,ratjbmx,
. stmx,psipol,chipol,real(index_rt) . stmx,psipol,chipol,real(index_rt),ajmxfi,dpdvmx
do i=1,nd do i=1,nd
if (ipec.eq.0) then if (ipec.eq.0) then