From 8c9273d684f19aa72e9fe71fde2d4437cec74881 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Wed, 28 Nov 2012 18:33:02 +0000 Subject: [PATCH] added reflection on limiter with arbitrary R,z shape --- Makefile | 4 +- src/gray.f | 197 +++++++++++++++++++++++++++++++++-------------------- 2 files changed, 126 insertions(+), 75 deletions(-) diff --git a/Makefile b/Makefile index 593ee12..173e322 100644 --- a/Makefile +++ b/Makefile @@ -2,7 +2,8 @@ EXE=gray # Objects list -OBJ=gray.o grayl.o green_func_p.o const_and_precisions.o itm_constants.o itm_types.o +OBJ=gray.o grayl.o reflections.o green_func_p.o \ + const_and_precisions.o itm_constants.o itm_types.o # Alternative search paths vpath %.f90 src @@ -21,6 +22,7 @@ $(EXE): $(OBJ) $(FC) $(FFLAGS) -o $@ $^ # Dependencies on modules +gray.o: green_func_p.o reflections.o green_func_p.o: const_and_precisions.o const_and_precisions.o: itm_types.o itm_constants.o itm_constants.o: itm_types.o diff --git a/src/gray.f b/src/gray.f index d2e3729..1d92780 100644 --- a/src/gray.f +++ b/src/gray.f @@ -225,7 +225,7 @@ c call gwork(j,k) c if(ierr.gt.0) then - print*,' IERR = ', ierr +! print*,' IERR = ', ierr if(ierr.eq.97) then c igrad=0 c ierr=0 @@ -257,11 +257,14 @@ c ierr=0 if(i.gt.1.and.psinv.ge.0.and.psinv.lt.psdbnd) . iov(j,k)=1 - if(iov(j,k).eq.1.and.psinv.ge.psdbnd) iov(j,k)=2 + if(iov(j,k).eq.1.and. + . (psinv.ge.psdbnd.or. + . (psinv.lt.1.0d0.and.(zzm.lt.zbmin.or.zzm.gt.zbmax)))) + . iov(j,k)=2 c iov=0 initially, iov=1 first entrance in plasma, c iov=2 first exit from plasma, iov=3 after 2nd entrance into plasma - if(index_rt.eq.1) then + if(index_rt.eq.1) then if(j.eq.1) then psinv11=psinv if(ipolc.eq.0.and.iov(j,k).eq.1) then @@ -278,10 +281,10 @@ c print*,' ' c print*,'Input coupled power (MW) =',p0mw c print*,' ' end if - if (ipass.gt.1) then - if(ipolc.eq.1.and.iov(j,k).eq.2.and.rrm.le.rcen) then + if (iov(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1 + . .and.ipolc.eq.1) then call pol_limit(qqout,uuout,vvout) - ipolc=2 + ipolc=2 call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl) strfl11=i*dst+dstvac call pol_limit(qqin2,uuin2,vvin2) @@ -300,7 +303,6 @@ c print*,' ' yyrfl(j,k,6)=anvrfl(3) tau1v(j,k)=tauv(j,k,iiv(j,k)) end if - end if else if(iov(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1) then call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl) @@ -537,6 +539,8 @@ c character*24 filenmeqq,filenmprf,filenmbm parameter(qe=4.8032d-10,me=9.1095d-28,vc=2.9979d+10) parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0) + parameter(nmx=8000,nbb=1000) + real*8 rlim(nbb),zlim(nbb) c common/xgcn/xgcn @@ -556,6 +560,7 @@ c common/igrad/igrad common/ipass/ipass common/rwallm/rwallm + common/limiter/rlim,zlim,nlim common/iieq/iequil common/icocos/icocos common/ixp/ixp @@ -827,6 +832,22 @@ c versus psi, rhop, rhot if (iequil.eq.1) call surf_anal + if (iequil.ne.2.or.ipass.lt.0) then +c set simple limiter as two cylindrical walls at rwallm and r00 + nlim=5 + rlim(1)=rwallm + rlim(2)=r00*1.d-2 + rlim(3)=rlim(2) + rlim(4)=rlim(1) + rlim(5)=rlim(1) + zlim(1)=-dst*nmx*1.d-2 + zlim(2)=zlim(1) + zlim(3)=dst*nmx*1.d-2 + zlim(4)=zlim(3) + zlim(5)=zlim(1) + ipass=abs(ipass) + end if + nfil=78 open(nfil,file='headers.txt',status= 'unknown') @@ -1047,6 +1068,7 @@ c dimension ffprim(nnw),pprim(nnw) dimension psi(nnw,nnh),rv(nnw),zv(nnh),psin(nnw,nnh),psinr(nnw) dimension rbbbs(nbb),zbbbs(nbb) + dimension rlim(nbb),zlim(nbb) c parameter(nrest=nnw+4,nzest=nnh+4) parameter(lwrk=4*(nnw+nnh)+11*(nrest+nzest)+nrest*nnh+nzest+54) @@ -1092,6 +1114,7 @@ c common/rhotsx/rhotsx common/rrtor/rrtor common/cnt/rup,zup,rlw,zlw + common/limiter/rlim,zlim,nlim c c read from file eqdsk c (see http://fusion.gat.com/efit/g_eqdsk.html) @@ -1267,41 +1290,45 @@ c call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz, . ffvpsi,cc11,lw11,iwrkd,ldiwrk,ier) c - if(ixp.ne.0) then - read (neqdsk,*) nbbbs,limitr - if(nbbbs.gt.0) then - if(ipsinorm.eq.1) - . read (neqdsk,*) (rbbbs(i),zbbbs(i),i=1,nbbbs) - if(ipsinorm.eq.0) - . read (neqdsk,2020) (rbbbs(i),zbbbs(i),i=1,nbbbs) - end if + read (neqdsk,*) nbbbs,nlim + if(nbbbs.gt.0) then + if(ipsinorm.eq.1) + . read (neqdsk,*) (rbbbs(i),zbbbs(i),i=1,nbbbs) + if(ipsinorm.eq.0) + . read (neqdsk,2020) (rbbbs(i),zbbbs(i),i=1,nbbbs) + end if + if(nlim.gt.0) then + if(ipsinorm.eq.1) + . read (neqdsk,*) (rlim(i),zlim(i),i=1,nlim) + if(ipsinorm.eq.0) + . read (neqdsk,2020) (rlim(i),zlim(i),i=1,nlim) + end if c c compute max and min z of last closed surface c + rbmin=rax + rbmax=rax + if (nbbbs.gt.1) then zbmin=1.0d+30 zbmax=-1.0d+30 - if (nbbbs.gt.1) then - do i=1,nbbbs - if(zbbbs(i).le.zbmin) then - zbmin=zbbbs(i) - rbmin=rbbbs(i) - end if - if(zbbbs(i).ge.zbmax) then - zbmax=zbbbs(i) - rbmax=rbbbs(i) - end if - end do - end if - if(zbmin.eq.zmnm) zbmin=zbmin+dz - if(rbmin.eq.rmnm) rbmin=rbmin+dr - if(zbmax.eq.zmxm) zbmax=zbmax-dz - if(rbmax.eq.rmxm) rbmax=rbmax-dr + do i=1,nbbbs + if(zbbbs(i).le.zbmin) then + zbmin=zbbbs(i) + rbmin=rbbbs(i) + end if + if(zbbbs(i).ge.zbmax) then + zbmax=zbbbs(i) + rbmax=rbbbs(i) + end if + end do else - zbmin=zmnm+dz - rbmin=rmnm+dr - zbmax=zmxm-dz - rbmax=rmxm-dr + zbmin=-1.0d+30 + zbmax=1.0d+30 end if + if(zbmin.le.zmnm) zbmin=zbmin+dz + if(rbmin.le.rmnm) rbmin=rbmin+dr + if(zbmax.ge.zmxm) zbmax=zbmax-dz + if(rbmax.ge.rmxm) rbmax=rbmax-dr c c scaling of f_poloidal c @@ -6332,9 +6359,9 @@ c computation of reflection coordinates and normal to the wall if(ivac.lt.0) then irfl=0 - xvrfl=xvout - xv=xvout - anvrfl=anv + xvrfl=xvout + xv=xvout + anvrfl=anv return end if @@ -6389,15 +6416,18 @@ c wave vector and electric field after reflection in lab frame end subroutine vacuum_rt(xvstart,xvend,ivac) + use reflections, only : inters_linewall,inside implicit none integer*4 ivac - real*8 x00,y00,z00,r00 - real*8 st,rs,rrm,psinv,rwallm,pi,dst,psdbnd,dstvac,deltawall + integer nbb,nlim,i,imax + parameter(nbb=1000) + real*8 st,rrm,zzm,psinv,dst,psdbnd,dstvac,smax real*8 anv(3),xvstart(3),xvend(3),walln(3),y(6),dery(6) + real*8 xv0(3) + real*8 rlim(nbb),zlim(nbb) - common/rwallm/rwallm common/wrefl/walln - common/mirr/x00,y00,z00 + common/limiter/rlim,zlim,nlim common/anv/anv common/dsds/dst common/psival/psinv @@ -6406,40 +6436,59 @@ c wave vector and electric field after reflection in lab frame c ivac=1 first interface plasma-vacuum c ivac=2 second interface vacuum-plasma after wall reflection c ivac=-1 second interface vacuum-plasma WITHOUT wall reflection + +! the real argument associated to xvstart is in a common block +! used by fwork!!! + xv0=xvstart + call inters_linewall(xv0/1.d2,anv,rlim(1:nlim),zlim(1:nlim), + . nlim,smax,walln) + smax=smax*1.d2 + rrm=1.d-2*sqrt(xv0(1)**2+xv0(2)**2) + zzm=1.d-2*xv0(3) + if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then + ! first wall interface is outside-inside + if (dot_product(walln,walln)