; FILE: bldyn_anim_1.idl [IDL source code] ;=========================================================== ; NONLINEAR AXISYMMETRIC BABCOCK-LEIGHTON DYNAMO ;----------------------------------------------------------- ; APAS7500 Fall 1997 Problem II.X.Y ; see Section X.Y of notes, and Appendix X ;----------------------------------------------------------- ;----------------------------------------------------------- ;----------------------------------------------------------- ; INPUT: this code requires an input file, the name of which ; being set immediately below ;=========================================================== pro animation ; solfile = 'saasfe/VarCo/Re=10Ca=10Co=8840Eta=100Ipar=2.i3e' ; solfile = 'saasfe/VarEta/Re=10Ca=10Co=17680Eta=1000Ipar=2.i3e' solfile = 'saasfe/Re=15Ca=35Co=5000Eta=100Ipar=0.i3e' titlelab ='!17a0 C!D!7a!n!17=+5 !7Dg=1/300 !17Rm=840' imakemovie = 0 device, DECOMPOSED = 0 parity = -1. makecolortable tvlct,rr,gg,bb,/get a=findgen(16)*(!PI*2./16.) usersym,cos(a),sin(a),/fill !ORDER = 1 imovie = 1 ;-------si animate=1 on va afficher xinteranimate----------- animate = 1 ;***** you shouldn't have to touch anything below this point ***** !ORDER = 1 close,1 & openr,1,/f77_unformatted,solfile title0=string(' ',format='(a80)') & ctime=1.e0 readu,1,title0 & titlelab=title0 print,'Using Input file: ',solfile print,'Run title: ',title0 idum=lonarr(6) vdum=fltarr(8) readu,1,idum numnp =idum(0) numel =idum(1) nelx =idum(2) nely =idum(3) nt0 =idum(4) nsteps=idum(5) readu,1,comega,calpha,reynolds,etadf,ybl,rzc,to title0 = "Standard Re="+strtrim(uint(reynolds),2)+" C!7a!6="+strtrim(fix(calpha),2)+" C!7X!6="+strtrim(long(comega),2)+" !7n!6df="+strtrim(uint(etadf),2) print,comega,calpha,reynolds,etadf,ybl,rzc,to nt0 = nt0*nsteps nnodes = numnp/2 nskip = 500 ; number of time steps skipped over; for sim6 et sim9 nt = 500 ; number fo time steps in animation xx = fltarr(nelx+1) yy = fltarr(nely+1) b_r = fltarr(nt+1,nelx+1) b_y = fltarr(nt+1,nelx+1) omg = fltarr(nelx+1,nely+1) domgdm = fltarr(nelx+1,nely+1) domgdr = fltarr(nelx+1,nely+1) eta = fltarr(nelx+1,nely+1) ur = fltarr(nelx+1,nely+1) uth = fltarr(nelx+1,nely+1) aquad = fltarr(nelx+1) daqdt = fltarr(nelx+1) bquad = fltarr(nelx+1,nely+1) b_y = fltarr(nelx+1,nely+1) to = to*nt/nt0 time = findgen(nt+1)/(nt)*to emag = fltarr(nt+1) aphiphi = fltarr(nelx+1,nely+1) aphitht = fltarr(nelx+1,nely+1) aphir = fltarr(nelx+1,nely+1) athtphi = fltarr(nelx+1,nely+1) athttht = fltarr(nelx+1,nely+1) athtr = fltarr(nelx+1,nely+1) arphi = fltarr(nelx+1,nely+1) artht = fltarr(nelx+1,nely+1) arr = fltarr(nelx+1,nely+1) ;---dérivé en x du tenseur alpha -- daphiphidx = fltarr(nelx+1,nely+1) daphithtdx = fltarr(nelx+1,nely+1) daphirdx = fltarr(nelx+1,nely+1) dathtphidx = fltarr(nelx+1,nely+1) dathtthtdx = fltarr(nelx+1,nely+1) dathtrdx = fltarr(nelx+1,nely+1) darphidx = fltarr(nelx+1,nely+1) darthtdx = fltarr(nelx+1,nely+1) darrdx = fltarr(nelx+1,nely+1) ;---dérivé en y du tenseur alpha -- daphiphidy = fltarr(nelx+1,nely+1) daphithtdy = fltarr(nelx+1,nely+1) daphirdy = fltarr(nelx+1,nely+1) dathtphidy = fltarr(nelx+1,nely+1) dathtthtdy = fltarr(nelx+1,nely+1) dathtrdy = fltarr(nelx+1,nely+1) darphidy = fltarr(nelx+1,nely+1) darthtdy = fltarr(nelx+1,nely+1) darrdy = fltarr(nelx+1,nely+1) ;---------------------------------- rotdiff = fltarr(nt+1,nelx-1) source = fltarr(nt+1,nelx-1) summ = fltarr(nt+1,nelx-1) readu,1,xx readu,1,yy readu,1,omg ; readu,1,domgdm ; readu,1,domgdr readu,1,ur readu,1,uth readu,1,eta ; ----------on lit alpha------------- readu,1,aphiphi ; readu,1,aphitht ; readu,1,aphr ; readu,1,athtphi ; readu,1,athttht ; readu,1,athtr ; readu,1,arphi ; readu,1,artht ; readu,1,arr ; ; ------on lit dérivé en x alpha----- ; readu,1,daphiphidx ; readu,1,daphithtdx ; readu,1,daphirdx ; readu,1,dathtphidx ; readu,1,dathtthtdx ; readu,1,dathtrdx ; readu,1,darphidx ; readu,1,darthtdx ; readu,1,darrdx ; ; ------on lit dérivé en y alpha----- ; readu,1,daphiphidy ; readu,1,daphithtdy ; readu,1,daphirdy ; readu,1,dathtphidy ; readu,1,dathtthtdy ; readu,1,dathtrdy ; readu,1,darphidy ; readu,1,darthtdy ; readu,1,darrdy xxa = acos(xx) yc = fltarr(nelx+1,nely+1) xc = fltarr(nelx+1,nely+1) ;--------- construct grid in real [x,y] space for i = 0,nely do begin for j = 0,nelx do begin sin1 = sqrt(1.-xx(j)^2) xc(j,i) = yy(i)*sin1 yc(j,i) = yy(i)*xx(j) endfor endfor ;---------- find boundaries of layer and top iblay = 0 ilay = 0 iint = 0 for k = 1,nely-1 do if yy(k+1) gt rzc and yy(k) le rzc then iint=k for k = 1,nely-1 do if yy(k+1) gt (rzc-ybl) and yy(k) le (rzc-ybl) then iblay=k for k = 1,nely-1 do if yy(k+1) gt 1.0 and yy(k) le 1.0 then itop=k print, iblay,yy(iblay) print, iint,yy(iint) print, itop,yy(itop) ;----------------------------------------------------------------------- ; select contour interval and color levels ; 1. choose contour interval and color levels cin = [ 0, 10, 11, 12, 13, 14, 15, 16, 17, $ 18, 19, 20, 21, 22, 23, 24, 25, $ 26, 27, 28, 29, 30 ] + 1 ; levb_y=[-10.,-3.e0,-1.0e-0,-3.e-1,-1.0e-1,-3.e-2,-1.0e-2,$ ; -3.e-3,-1.e-3,-3.e-4,-1.e-4, 1.e-4, 3.e-4, 1.0e-3, 3.e-3,$ ; 1.0e-2, 3.e-2, 1.0e-1, 3.e-1, 1.0e0, 3.e0, 10.] ;levb_y=[-100.,-20.,-9.,-8.,-7.,-6.,-5.,-4.,$ ; -3.,-2.,-1., 1., 2., 3., 4.,$ ; 5., 6., 7., 8., 9., 10., 100.] levb_y = [-1.1,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,$ 0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1] ; levb_y = [-0.9,-0.85,-0.8,-0.75,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,$ ; 0.2,0.3,0.4,0.5,0.6,0.7,0.75,0.8,0.85,0.9]/10. levquad = levb_y ;----------------------------------------------------------------------------- ; prepare output file for movie if imovie eq 1 then begin ; workstation animation nnx = 350 nny = 350 xr = [0.,1.1] yr = [0.,1.1] if xx[0] lt -0.5 then begin yr = [-1.,1.] nny = 700 endif pos0 = [0.1,0.05,0.975,0.95] mimage = bytarr(nnx,nny,nt+1) itr = 0 xlab = 0.3 ylab = -0.05 xlab2 = 0. ylab2 = -0.1 thk = 2 chsiz = 1.5 ncyc = 1. ccol = [7,7,7,7,7,7,7,7,7,7,7,3,3,3,3,3,3,3,3,3,3,3] ccol2 = [250,123,12,35,67,53,23,12,23,44] lev2 = findgen(200)/10 - 10 endif window,1,retain=2,xpos=50,ypos=200,xs=nnx,ys=nny ; Initialize XINTERANIMATE: if animate eq 1 then xinteranimate, set=[nnx, nny, nt+1], /mpeg_open , mpeg_bitrate= 104857200 relab = '!6Re(s)' imlab = '!6Im(s)' colab = '!6C!7X!6' calab = '!6C!7a!6' colab2 = '!6=1.E5' calab2 = '!6='+string(calpha,format='(f4.1)') ;----------------------------------------------------------------------------- ; loop over timesteps vdum = fltarr(5) ifr = -1 ;---------- read results for eigenvector and eigenvalue ; (skip over nskip solutions)window for iskip = 0, nskip-1 do begin readu,1,vdum readu,1,bquad readu,1,b_y endfor for it = 0, nt do begin readu,1,vdum time(it) = vdum(0) & emag(it)=vdum(2) readu,1,bquad readu,1,b_y ; normalization if it eq 0 then begin aqmax = max(abs(bquad)) bymax = max(abs(b_y)) endif bquad = bquad/aqmax b_y = b_y/bymax ;------------------------------------------------------------------------------ !X.STYLE = 5 !Y.STYLE = 5 ;---------- 1. plot color-filled contours for B-phi ; b_y *= 20. ; bquad *= 20. if ifr eq 0 then erase contour,pos = pos0,levels=levb_y,/fill,c_color=cin,xrange=xr,yrange=yr,$ b_y(itr:nelx,0:itop),xc(itr:nelx,0:itop), yc(itr:nelx,0:itop) oplot,/noclip,xc(itr:nelx,itop:itop), yc(itr:nelx,itop:itop),thick=thk oplot,/noclip,xc(itr:nelx,iint:iint), yc(itr:nelx,iint:iint),thick=2,$ linestyle = 2 if xx[0] lt -0.5 then oplot,/noclip,[0.,0.],[0.,-1.],thick=thk oplot,/noclip,[0.,0.],[0.,1.],thick=thk oplot,/noclip,[0.,1.],[0.,0.],thick=thk ;---------- 2. overplot contours for B-pol contour,/follow,/noclip,/noerase,pos=pos0,levels=levquad,$ xrange = xr,yrange=yr,$ c_labels = [0,0,0,0,0,0,0,0,0],c_color=ccol,c_thick=thk,$ bquad(itr:nelx,0:itop+8),xc(itr:nelx,0:itop+8), yc(itr:nelx,0:itop+8) timelab = '!8t/!7s!17='+string(time(it),format='(f6.4)') xyouts,xlab2,ylab ,charsize=chsiz,timelab,color=5 xyouts,charsize=1.5,0.,1.05,title0 ; 5. color scale colorscale,pos0,levb_y,cin,-0.1,-0.05,0.,1.0,thk ; 6. store imagesolfile='saasfe/test2_corinne.i3e' ifr = ifr+1 mimage2 = bytarr(nnx,nny,3,nt+1) if imovie eq 1 then begin mimage2(*,*,*,ifr:ifr) = tvrd(0,0,nnx,nny,1,true=2) mimage(*,*,ifr:ifr) = tvrd(0,0,nnx,nny,1) endif ; --- xinteranimate ------------- ; Load the images into XINTERANIMATE: if animate eq 1 then xinteranimate, frame = it, window = 1 ; -- store animation as animated gif if imakemovie eq 2 then begin if it eq 0 then mpegID = MPEG_OPEN([nnx,nny]) colordata = tvrd() image24 = bytarr(3,nnx,nny) tvlct,r,g,b,/get image24[0,*,*] = r(colordata[*,*]) image24[1,*,*] = g(colordata[*,*]) image24[2,*,*] = b(colordata[*,*]) MPEG_PUT,mpegID,image=image24,FRAME=it endif ; -- store animation as mpeg file if imakemovie eq 1 then begin if it eq 0 then mpegID = MPEG_OPEN([nnx,nny]) colordata = tvrd() image24 = bytarr(3,nnx,nny) tvlct,r,g,b,/get image24[0,*,*] = r(colordata[*,*]) image24[1,*,*] = g(colordata[*,*]) image24[2,*,*] = b(colordata[*,*]) MPEG_PUT,mpegID,image=image24,FRAME=it endif endfor if imakemovie eq 1 then begin MPEG_SAVE,mpegID,FILENAME='animBLgif.mpg' MPEG_CLOSE,mpegID endif ; Play the animation: if animate eq 1 then xinteranimate, /keep_pixmaps, /MPEG_CLOSE print,max(b_y),max(bquad) end