From 1a60867b5b1679d36f02c337eb37a998f87a2a0d Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Mon, 27 Jan 2014 15:38:17 +0000 Subject: [PATCH] pec subroutine: new algorithm to compute Jphi --- scripts/beta-alpha_beam.sh | 12 +-- src/gray.f | 152 ++++++++++++++++++++++++++++++++----- 2 files changed, 141 insertions(+), 23 deletions(-) diff --git a/scripts/beta-alpha_beam.sh b/scripts/beta-alpha_beam.sh index cf32f71..76a3292 100755 --- a/scripts/beta-alpha_beam.sh +++ b/scripts/beta-alpha_beam.sh @@ -1,9 +1,11 @@ #!/bin/bash # set input/working/output folders -in=/Users/daniela/Desktop/scenario_001/files5s -work=/Users/daniela/Desktop/scenario_001/tmpab_beam -out=/Users/daniela/Desktop/scenario_001/outab_beam_550 +in=/Users/lorenzo/ITER/UL-Grant-161/performance-analysis/case006/in +work=/Users/lorenzo/Desktop/tmp/sc6 +out=/Users/lorenzo/Public/case006-mappe-ab-beam + +scen='006' # launching mirrors mirror=( "USM" "LSM" ) @@ -58,8 +60,8 @@ cp $in/${base}.prf ${base}.prf #loop on USM/LSM mirrors for i in 0 1; do # set output file names - id="${mirror[i]}_t$t" - f7file=$out/f7_${id}.txt + id="${mirror[i]}_Scen${scen}_t$t" + f7file=$out/ECCD_vs_ba_${id}.txt f48file=$out/f48_${id}.txt rm -f $out/log_${id}.txt # add header in output file diff --git a/src/gray.f b/src/gray.f index 6f368da..90b3aac 100644 --- a/src/gray.f +++ b/src/gray.f @@ -512,9 +512,9 @@ c 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(12,*) ' #i sst psi w1 w2' - write(7,*)'#Icd Pa Jphimx dPdVmx '// - .'rhotj rhotjava rhotp rhotpav '// - .'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt' + write(7,*)'#Icd Pa Jphip dPdVp '// + .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav '// + .'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx' write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt' else @@ -1118,6 +1118,7 @@ c common/cofffp/cfp common/fpas/fpolas common/rhotsx/rhotsx + common/phitedge/phitedge common/rrtor/rrtor common/cnt/rup,zup,rlw,zlw common/limiter/rlim,zlim,nlim @@ -1457,8 +1458,10 @@ c compute B_toroidal on axis c compute normalized rho_tor from eqdsk q profile call rhotor(nr) -c phitedge=abs(psia)*rhotsx*2*pi -c rrtor=sqrt(phitedge/abs(btrcen)/pi) + phitedge=abs(psia)*rhotsx*2*pi + rrtor=sqrt(phitedge/abs(btrcen)/pi) + call rhopol +c print*,rhotsx,phitedge,rrtor,abs(psia) c compute flux surface averaged quantities @@ -1822,6 +1825,7 @@ c common/qpsi/qpsi common/rhotsx/rhotsx common/crhot/crhot + common/cq/cq c c normalized toroidal rho : ~ Integral q(psi) dpsi c @@ -1847,6 +1851,21 @@ c return 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) implicit real*8(a-h,o-z) parameter(nnw=501) @@ -1883,7 +1902,57 @@ c if(ip.eq.nintp) ip=nintp-1 dps=rpsi-rpstab(ip) 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 end @@ -2158,6 +2227,8 @@ c dimension pstab(nnintp),varea(nnintp),vvol(nnintp),rpstab(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 vcurrp(nnintp),vajphiav(nnintp),qqv(nnintp) dimension rcon(2*ncnt+1),zcon(2*ncnt+1) @@ -2186,6 +2257,9 @@ c common/cnt/rup,zup,rlw,zlw common/bound/zbmin,zbmax common/rarea/rarea + common/phitedge/phitedge + common/cdadrhot/cdadrhot + common/cdvdrhot/cdvdrhot c common/coffvl/ch common/coffdvl/ch01 @@ -2368,6 +2442,11 @@ c ratio_pltor = Jcd_||/J_phi Jcd_|| = vratjb(jp)=ratio_cdbtor qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi)) 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 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/ end do write(56,*)' #psi rhot_eq rhot_av || |Bmx| |Bmn|'// - .' Area Vol |I_pl| qq fc ratio_cdbtor' + .' Area Vol |I_pl| qq fc' qqv(1)=qqv(2) vajphiav(1)=vajphiav(2) @@ -2467,6 +2546,10 @@ c spline coefficients of rhot call difcs(rpstab,ffc,nintp,iopt,cfc,ier) iopt=0 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 @@ -2485,9 +2568,37 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda return 99 format(20(1x,e12.5)) end -c -c -c + + function fdadrhot(rpsi) + 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 implicit real*8 (a-h,o-z) parameter(jmx=31,kmx=36,nmx=8000) @@ -6180,7 +6291,12 @@ c of gaussian profile ajcdbv(i)=ajphiv(i)*ratjbv(i) ajplv(i)=ajphiv(i)*ratjplv(i) 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, . drhotp,drhopp) if(ieccd.ne.0) then @@ -6253,18 +6369,18 @@ c dPdV [MW/m^3], Jcd [MA/m^2] currtka=currt*1.0d3 write(6,*)' ' - write(6,*)'#beta0 alpha0 Icd Pa Jphimx dPdVmx '// - .'rhotj rhotjava rhotp rhotpav '// - .'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt' - write(6,99) beta0,alpha0,currtka,pabstot,ajmxfi,dpdvmx, + write(6,*)'#beta0 alpha0 Icd Pa Jphip dPdVp '// + .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav '// + .'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx' + write(6,99) beta0,alpha0,currtka,pabstot,ajphip,dpdvp, . rhotjfi,rhotjava,rhotp,rhotpav, . 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, . drhotjava,drhotpav,ratjbmx, - . stmx,psipol,chipol,real(index_rt) + . stmx,psipol,chipol,real(index_rt),ajmxfi,dpdvmx do i=1,nd if (ipec.eq.0) then