program dyn2d c********************************************************************** c Axisymmetric Mean-field Dynamo in Spherical Geometry c APAS7500, Fall 97, Problems IV.5.1 through IV.5.5 c Written by Paul Charbonneau, March 97 c======================================================================= c This codes solves the dynamo eigenvalue problem for differential c rotation and alpha-effect profiles provided by the user, in the c form of two subroutines (ALPXY, OMGXY, see below). c c The code works on the scaled version of the equations c (see \S IV.4.2 of the class notes) and matches the solution to c a potential field in the exterior (r/R>1), assuming constant c magnetic diffusivity in the interior (r/R<1). c c Spatial discretization makes use of bilinear finite elements, c on a spatial mesh of fixed size N_r=N_theta=64x48 c The eigenvalue problem is solved by inverse iteration. c====================================================================== c OUTPUT: This code produces a single output file called c c solution.i3e c c This output file is IEEE Fortran unformatted; it is meant to be read c by the three IDL subroutines c c animsol.idl <-- animation of the solution over one cycle c medplane.idl <-- plots solution for six phases over one cycle c butterfly.idl <-- plots butterfly diagrams for B_phi and B_r c====================================================================== c INPUT REQUIRED FROM USER: c c two fortran functions; examples are provided at the very end c of this file c c ALPXY --> computes alpha effect functional at (cos(theta),r/R) c OMGXY --> computes angular velocity at (cos(theta),r/R) c c Four scalar quantities, values to be set at beginning of c executable section of the main program (15 lines down from here!) c c comega --> is the C_Omega dynamo number c calpha --> is the C_Omega dynamo number c reeig0 --> is the real part of the trial eigenvalue c imeig0 --> is the imaginary part of the trial eigenvalue c c If you are doing problem IV.5.4, you will also need to alter c the fortran function c c DMAGXY --> computes magnetic diffusivity at (cos(theta),r/R) c c********************************************************************** real imeig0 complex teigvl common/probt/ title character*64 title c------------------------------------------------------------------------ c --> unit 8: output open(8,file='solution.i3e',form='unformatted') c---------- define problem data title='Test case, from Stix 1976 [Model 2] alpha domg/dr>0' c you MUST set values for the following: c --> comega is the C_Omega dynamo number c --> calpha is the C_Omega dynamo number c --> reeig0 is the real part of the trial eigenvalue c --> imeig0 is the imaginary part of the trial eigenvalue c the following sets a flag that enforce a parity on the c solutions. DO NOT TOUCH unless you are doing problem IV.5.5 ipar=0 c following for omg1, alp1, C_alpha<0 calpha=-120. comega=120. reeig0=0. imeig0=100. c=============== LEAVE AS IS BEYOND THIS POINT ====================== teigvl=cmplx(reeig0,imeig0) c---------- generate mesh call mesh2d(ipar) c---------- set up initial guess for eigen- value/vector call eigvv0(calpha,comega,teigvl) c---------- assemble system matrices call form c---------- apply boundary conditions call aplybc c---------- solver call calfun c---------- output call writout c==================================================================== stop ' **END**' end c********************************************************************** c YOU SHOULD NOT HAVE TO TOUCH ANYTHING BEYOND THIS POINT * c********************************************************************** subroutine eigvv0(calph,comeg,teig) c sets up some parameters relating to inverse iteration parameter(nm=2145,ne=2048,nb=141,ng=4) complex teig common/scale/ calpha,comega common/nodes/ x(nm),y(nm),numnp common/elmnts/ ntype(ne),ngpint(ne),icon(ne,4),numel common/syseqs/ ka(nb,2*nm),kb(nb,2*nm),kab(nb,2*nm), + fb(2*nm),q(2*nm),bca(2*nm) complex ka,kb,kab,fb,q,bca common/initq/ qsav(2*nm) complex qsav common/eigval/ neig,teigvl,eigvl,qnm1(2*nm) complex eigvl,qnm1,teigvl common/probt/ title character*64 title data pi/3.1415926536/ c---------------------------------------------------------------------- calpha=calph comega=comeg teigvl=teig c---------- 2. guess at eigenvector: c use eigenvector from previous sequence do 10 i=1,numnp qnm1(i)=qsav(i) 10 continue write(*,1) title,calph,comeg,real(teigvl),aimag(teigvl) return c------------------------------------------------------------------------- 1 format(' ',/,64('*'),/, + '* LINEAR SPHERICAL MEAN-FIELD DYNAMO',/,64('*'),/, + 1x,a64,//, + 1x,'Input Parameters:',/, + 1x,'C_alpha =',e12.4,/, + 1x,'C_omega =',e12.4,/, + 1x,'Re(Trial eigenvalue) =',e12.4,/, + 1x,'Im(Trial eigenvalue) =',e12.4,/) end c********************************************************************** subroutine writout parameter(nm=2145,ne=2048,nb=141,ng=4) common/syseqs/ ka(nb,2*nm),kb(nb,2*nm),kab(nb,2*nm), + fb(2*nm),q(2*nm),bca(2*nm) complex ka,kb,kab,fb,q,bca common/nodes/ x(nm),y(nm),numnp common/elmnts/ ntype(ne),ngpint(ne),icon(ne,4),numel common/gnode/ yy(nm-ne),xx(nm-ne),nely,nelx,nnodes common/scale/ calpha,comega common/eigval/ neig,teigvl,eigvl,qnm1(2*nm) complex eigvl,qnm1,teigvl common/iter/ itermax,icnvrg,toler,rayquo dimension qi1(nm),qr1(nm),qi2(nm),qr2(nm) common/probt/ title character*64 title data pi/3.1415926536/ c------------------------------------------------------------------------ write(8) title write(8) numnp,numel,nelx,nely write(8) calpha,comega write(8) (xx(j),j=1,nelx+1) write(8) (yy(j),j=1,nely+1) c--------- separate real and imaginary parts, for both displacements do 1 i=1,numnp/2 qr1(i)= real(q(2*i-1)) qi1(i)=aimag(q(2*i-1)) qr2(i)= real(q(2*i)) qi2(i)=aimag(q(2*i)) 1 continue write(8) eigvl,rayquo write(8) (qr1(j),j=1,numnp/2) write(8) (qi1(j),j=1,numnp/2) write(8) (qr2(j),j=1,numnp/2) write(8) (qi2(j),j=1,numnp/2) return c-------------------------------------------------------------------------- end c************************************************************************** subroutine initdt c reads initial data, sets limits and integration points parameter(nm=2145,ne=2048,nb=141,ng=4) common/elspec/ nnodes(30),intdef(30),ioptau(30) common/glint/ gp(5,5),wt(5,5) c------------------------------------- do 10 i=1,30 nnodes(i)=0 intdef(i)=0 ioptau(i)=0 10 continue nnodes(21)=3 nnodes(22)=6 nnodes(23)=4 nnodes(24)=8 nnodes(25)=4 nnodes(26)=4 intdef(21)=0 intdef(22)=4 intdef(23)=2 intdef(24)=3 intdef(25)=2 intdef(26)=2 ioptau(21)=1 ioptau(22)=4 ioptau(23)=1 ioptau(24)=4 ioptau(25)=1 ioptau(26)=1 c---------- gauss points and weights for gauss-legendre quadrature do 50 n=1,5 do 50 l=1,5 gp(n,l)=0. wt(n,l)=0. 50 continue c---------- for one gauss point gp(1,1)=0. wt(1,1)=2. c---------- for two gauss points gp(2,1)=-1./sqrt(3.) gp(2,2)=-gp(2,1) wt(2,1)=1. wt(2,2)=1. c---------- for three gauss points gp(3,1)=-sqrt(3./5.) gp(3,2)=0. gp(3,3)=-gp(3,1) wt(3,1)=5./9. wt(3,2)=8./9. wt(3,3)=wt(3,1) c---------- for four gauss points gp(4,1)=-sqrt((15.+2.*sqrt(30.))/35.) gp(4,2)=-sqrt((15.-2.*sqrt(30.))/35.) gp(4,3)=-gp(4,2) gp(4,4)=-gp(4,1) wt(4,1)=49./(6.*(18.+sqrt(30.))) wt(4,2)=49./(6.*(18.-sqrt(30.))) wt(4,3)=wt(4,2) wt(4,4)=wt(4,1) 500 continue c----------------------------------------------------------------------- end c*********************************************************************** subroutine aplybc c applies the natural and essential boundary conditions to c the system matrices parameter(nm=2145,ne=2048,nb=141,ng=4) common/genrl/ mband,mfbw common/essbc/ uebc(2*ne),iuebc(2*ne),numebc complex uebc common/nodes/ x(nm),y(nm),numnp common/elmnts/ ntype(ne),ngpint(ne),icon(ne,4),numel common/syseqs/ ka(nb,2*nm),kb(nb,2*nm),kab(nb,2*nm), + fb(2*nm),q(2*nm),bca(2*nm) complex ka,kb,kab,fb,q,bca common/elvar/ xi,eta,jac(ne,ng),jacinv(ne,ng,2,2), + xx(4),yy(4),phi(4), + dphdx(4),dphdy(4),dphdxi(4),dphdet(4),ielno,nnnd real jacinv,jac complex czero,cone c******************** apply mixed boundary conditions, by modifying c diagonal of stiffness matrix c write(*,*) ' ***** entering aplybc *****' c******************** apply the essential boundary conditions, if any, c by imposing constraint equations nq=2 400 if(numebc.le.0) return czero=cmplx(0.,0.) cone=cmplx(1.,0.) c********** enforce dirichlet B.C. by applying constraint equations to system c---------- 1. modify system matrix do 470 ibc=1,numebc jj=iuebc(ibc) c write(*,*) ' APLYBC: working on condition #',ibc,' node #',jj c --> zero matrix "line" (shifted version) do 460 i=2,mband ishp=mband-i+1 ishm=mband+i-1 jrp=jj+i-1 jrm=jj-i+1 if(jrm.ge.1) bca(jrm)=bca(jrm)-kab(ishp,jj)*uebc(ibc) if(jrm.ge.1) kab(ishm,jrm)=czero if(jrp.le.numnp) bca(jrp)=bca(jrp)-kab(ishm,jj)*uebc(ibc) if(jrp.le.numnp) kab(ishp,jrp)=czero 460 continue c --> zero matrix column do 465 i=1,mfbw kab(i,jj)=czero 465 continue c --> set corresponding diagonal element of system matrix to Real(one) kab(mband,jj)=cone 470 continue c write(*,*) ' ***** exiting aplybc *****' return c------------------------------------ end c*********************************************************************** subroutine calfun c----------------------------------------------------------------------- c solves the eigenvalue problem using inverse iteration | c (see sects. 7.6.1 and 7.7.8 of Golub and Van Loan) | c | c Note kab(i,j) = [A] - mu[B] | c ka (i,j) = [A] | c kb (i,j) = [B] | c----------------------------------------------------------------------- parameter(nm=2145,ne=2048,nb=141,ng=4) common/genrl/ mband,mfbw common/nodes/ x(nm),y(nm),numnp common/elmnts/ ntype(ne),ngpint(ne),icon(ne,4),numel common/eigval/ neig,teigvl,eigvl,qnm1(2*nm) complex eigvl,qnm1,teigvl common/iter/ itermax,icnvrg,toler,rayquo common/essbc/ uebc(2*ne),iuebc(2*ne),numebc complex uebc common/syseqs/ ka(nb,2*nm),kb(nb,2*nm),kab(nb,2*nm), + fb(2*nm),q(2*nm),bca(2*nm) complex ka,kb,kab,fb,q,bca dimension qn(2*nm),fa(2*nm) complex qn,fa,emu,eign,eignm1,cka,ckb,suma,sumb,czero equivalence (qn,q) c----------------------------------- c******************** eigenvalue problem by inverse iteration c write(*,*) ' ***** entering calfun *****' c write(*,*) 'numnp = ',numnp c---------- save load vector and function at start of interval itermax=30 toler=1.e-6 nq=2 czero=cmplx(0.,0.) emu=teigvl eignm1=emu write(*,1) 1 format(1x,'BEGINNING INVERSE ITERATION',//, + 1x,'iter',2x,'Re(eig)',8x,'Im(eig)',8x,'Ray. Quo.') c---------- triangularize system matrix ( [A]-mu[B] ) call tribu(numnp,mband,mfbw,kab) c---------- iterate until tolerance is reached do 480 n=1,itermax c---------- compute [B] q(n-1) do 430 j=1,numnp cka=czero jpmb=j+mband do 431 i=max(1,jpmb-numnp),min(mfbw,jpmb-1) jshift=jpmb-i cka=cka+kb(i,jshift)*qnm1(jshift) 431 continue fb(j)=cka + bca(j) 430 continue c---------- solve ( [A]-mu[B] )q(n) = [B]q(n-1) for q(n) call rhsbu(numnp,mband,mfbw,kab,fb,qn) c---------- compute 2-norm of q(n), and normalize q(n) sum2n=0. do 440 j=1,numnp sum2n=sum2n + conjg(qn(j))*qn(j) 440 continue qn2n = sqrt(sum2n) c write(*,*) ' ***** 2-norm = ',qn2n qn2ni=1./qn2n do 441 j=1,numnp qn(j)=qn(j)*qn2ni 441 continue c---------- compute [A] q(n) and [B] q(n) do 470 j=1,numnp cka=czero ckb=czero jpmb=j+mband do 471 i=max(1,jpmb-numnp),min(mfbw,jpmb-1) jshift=jpmb-i cka=cka+ka(i,jshift)*qn(jshift) ckb=ckb+kb(i,jshift)*qn(jshift) 471 continue fa(j)=cka fb(j)=ckb 470 continue c H H c---------- (b) compute q(n) [A] q(n) and q(n) [B] q(n) suma=czero sumb=czero do 472 j=1,numnp suma=suma + conjg(qn(j))*fa(j) sumb=sumb + conjg(qn(j))*fb(j) 472 continue c---------- compute tentative eigenvalue eign = suma/sumb c---------- Check if tolerance is reached ennorm=conjg(eign)*eign rayquo=abs( (ennorm-conjg(eignm1)*eignm1)/ennorm ) write(*,902) n,eign,rayquo 902 format(i4,1p3e15.6) if(rayquo.le.toler) goto 500 c---------- save tentative eigen- value and vector for next iteration eignm1 = eign do 460 i=1,numnp qnm1(i)=qn(i) 460 continue c---------- go to next iteration 480 continue write(*,491) eign,eignm1 icnvrg=1 stop ' *FAIL*' 500 continue write(*,490) eign,eignm1 icnvrg=0 eigvl=eign return c------------------------------------- 490 format(//,' ITERATION FOR EIGENVALUE/VECTOR HAS CONVERGED', + //,5x,'eig(n-1) = ',2e16.8, + /,5x,'eig(n) = ',2e16.8) 491 format(//,' ITERATION FOR EIGENVALUE/VECTOR HAS *NOT* CONVERGED', + //,5x,'eig(n-1) = ',2e16.8, + /,5x,'eig(n) = ',2e16.8) end c****************************************** subroutine tribu(n,mhb,mb,gk) c****************************************** c triangularized shifted matrix c (modified alg., after Golub & Van Loan, p. 152) c****************************************** parameter(nm=2145,ne=2048,nb=141,ng=4) dimension gk(nb,2*nm) complex gk,c,cone,czero cone=cmplx(1.,0.) czero=cmplx(0.,0.) c write(*,*) ' ***** entering tribu *****' c write(*,*) n,mhb,mb do 1 k=1,n-1 c if(gk(mhb,k).eq.0.)goto 99 c if(gk(mhb,k).eq.1.)goto 5 c=cone/gk(mhb,k) do 2 i=mhb+1,mb gk(i,k)=c*gk(i,k) 2 continue 5 continue do 3 j=1,min(mhb-1,n-k) kpj=k+j c=gk(mhb-j,kpj) if(c.ne.czero)then CDIR$ IVDEP do 4 i=mhb+1,min(mb,n-k+mhb) imj=i-j gk(imj,kpj)=gk(imj,kpj) + -c*gk(i,k) 4 continue end if 3 continue 1 continue return 99 write(*,100) k,gk(mhb,k) stop 'tribu' 100 format(//,' set of equations may be singular..',/, +' diagonal term of equation',i5,' at tribu is equal to', +1p2e15.8) end c****************************************** subroutine rhsbu(n,mhb,mb,gk,gf,u) c****************************************** c forward reduction & backsubstitution c (modified alg., after Golub & Van Loan, p. 152) c****************************************** parameter(nm=2145,ne=2048,nb=141,ng=4) dimension gk(nb,2*nm),gf(2*nm),u(2*nm) complex gk,gf,u,cone,czero czero=cmplx(0.,0.) c---------- forward reduction do 1 j=1,n c if(gk(mhb,j).eq.czero)goto 99 mhbmj=mhb-j do 2 i=j+1,min(n,j+mhb-1) ishift=i+mhbmj gf(i)=gf(i)-gk(ishift,j)*gf(j) 2 continue 1 continue c---------- backsubstitution do 3 j=n,1,-1 c gf(j)=gf(j)/gk(mhb,j) u(j)=gf(j)/gk(mhb,j) mhbmj=mhb-j do 4 i=max(1,j-mhb+1),j-1 ishift=i+mhbmj c gf(i)=gf(i)-gk(ishift,j)*gf(j) gf(i)=gf(i)-gk(ishift,j)*u(j) 4 continue c u(j)=gf(j) 3 continue c write(*,*) ' counters: ic1, ic2 = ',ic1,ic2 return 99 write(*,100) j,gk(mhb,j) stop 'rhsbu' 100 format(//,' set of equations are singular.',/, +' diagonal term of equation',i5,' at tribu is equal to', +1p2e15.8) end c***************************************************************************** subroutine form c*************************************************************************** c form the system matrices from element matrices contributions * c * c axisymmetric problem --> TWO coupled equations * c*************************************************************************** parameter(nm=2145,ne=2048,nb=141,ng=4) common/genrl/ mband,mfbw common/syseqs/ ka(nb,2*nm),kb(nb,2*nm),kab(nb,2*nm), + fb(2*nm),q(2*nm),bca(2*nm) complex ka,kb,kab,fb,q,bca common/nodes/ x(nm),y(nm),numnp common/elmnts/ ntype(ne),ngpint(ne),icon(ne,4),numel common/eigval/ neig,teigvl,eigvl,qnm1(2*nm) complex eigvl,qnm1,teigvl common/elspec/ nnodes(30),intdef(30),ioptau(30) common/elstf/ aer(8,8),aei(8,8),ber(8,8),bei(8,8) common/elvar/ xi,eta,jac(ne,ng),jacinv(ne,ng,2,2), + xx(4),yy(4),phi(4), + dphdx(4),dphdy(4),dphdxi(4),dphdet(4),ielno,nnd real jacinv,jac complex czero c----------------------------------- c write(*,*) ' ***** entering form' c---------- initialize arrays for the system matrices nq=2 nqm1=nq-1 czero=cmplx(0.,0.) do 60 i=1,mfbw do 61 j=1,numnp kab(i,j)=czero ka(i,j) =czero kb(i,j) =czero bca(j) =czero 61 continue 60 continue c******************** loop over all elements nnd=4 do 500 ielno=1,numel c write(*,*) ' FORM: computing element matrix for element #',ielno c---------- form element matrices call elem23c c write(*,*) ' FORM: assembling element matrix for element #',ielno c---------- assemble the system matrices from the element matrices do 440 i=1,nnd i1 =nq*i-nqm1 ii1 =nq*icon(ielno,i)-nqm1 ii1pmb=ii1+mband do 419 j=0,nqm1 ij=ii1+j 419 continue do 420 j=1,nnd j1 =nq*j-nqm1 jj1=nq*icon(ielno,j)-nqm1 ish=ii1pmb-jj1 c --> alternate assembly : loop procedure do 421 k1=0,nqm1 ist=ish-k1 jk=jj1+k1 jae=j1+k1 do 422 k2=0,nqm1 iae=i1+k2 kab(ist+k2,jk) =kab(ist+k2,jk) + +cmplx(aer(iae,jae),aei(iae,jae)) + -teigvl*cmplx(ber(iae,jae),bei(iae,jae)) ka (ist+k2,jk) =ka (ist+k2,jk) + +cmplx(aer(iae,jae),aei(iae,jae)) kb (ist+k2,jk) =kb (ist+k2,jk) + +cmplx(ber(iae,jae),bei(iae,jae)) 422 continue 421 continue 420 continue 440 continue 500 continue c write(*,*) ' ***** exiting form' return end c****************************************************************************** subroutine header(ityp,ngp,incr,incw) parameter(nm=2145,ne=2048,nb=141,ng=4) common/gnode/ yy(nm-ne),xx(nm-ne),nely,nelx,nnodes common/mesht/ nprint,nt,xtime,times(500),nttot,nstp, + thetas(500),omegazc(500),domgdt(500) common/probt/ title character*64 title character*10 date,clock write(*,1) c write(*,2) date(),clock() write(*,3) title write(*,'()') if(incr.eq.0) write(*,*) ' Initial Run' if(incr.eq.1) write(*,*) ' Running from Restart File' if(incw.eq.0) write(*,*) ' Final Run: no Restart File' if(incw.eq.1) write(*,*) ' Creating Restart File' return c.............................................................................. 1 format(1x,60('*'),/,1x,'*',10x, +'KINEMATIC INTERFACE DYNAMOS',/,1x,'*',10x, +'KINEMATIC INDUCTION EQUATION SOLVER', +/,1x,60('*'),/) 2 format(1x,'submitted on',a10, + /,1x,' at',a10) 3 format(1x,'Run: ',a64,/) 4 format(1x,'Input Parameters:',/,5x, + i10, 5x,'ybl ',10x,'BL thickness',/,5x, + f10.2,5x,'eta2 ',10x,'Mag Reynolds number (>rzc)',/,5x, + f10.2,5x,'alpl ',10x,'Helicity length scale',/,5x, + f10.2,5x,'omega0 ',10x,'Initial Angular velocity',/,5x, + i10 ,5x,'nelx ',10x,'Number of elements in angular', + ' direction',/,5x, + i10 ,5x,'nelx ',10x,'Number of elements in radial', + ' direction') 5 format(5x,f10.2,5x,'concr ',10x,'Element concentration factor', + ' for radial direction',/,5x, + i10 ,5x,'nt ',10x,'Number of time step groups',/,5x, + e10.2,5x,'to ',10x,'Integration time (in yrs)',/,5x, + i10, 5x,'ityp ',10x,'Element type',/,5x, + i10, 5x,'ngp ',10x,'Number of points for gaussian', + ' quadrature',//) end c***************************************************************************** subroutine elem23c c***************************************************************************** c master routine for element 23c, a quadrilateral bilinear isoparametric * c element for a system of two coupled parabolic PDE * c***************************************************************************** parameter(nm=2145,ne=2048,nb=141,ng=4) common/elstf/ aer(8,8),aei(8,8),ber(8,8),bei(8,8) common/scale/ calpha,comega common/nodes/ x(nm),y(nm),numnp common/elmnts/ ntype(ne),ngpint(ne),icon(ne,4),numel common/intld/ fintu(ne),fintb(ne) common/elvar/ xi,eta,jac(ne,ng),jacinv(ne,ng,4), + xx(4),yy(4),phi(4), + dphdx(4),dphdy(4),dphdxi(4),dphdet(4),ielno,nnd real jacinv,jac common/phys/ xgp(ne,ng),ygp(ne,ng) common/glint/ gp(5,5),wt(5,5) common/gnode/ ydum(nm-ne),xdum(nm-ne),nely,nelx,indum real mu data pi/3.1415926536/ c********** form elements matrices and vectors c---------- initialisations do 20 i=1,4 jj=icon(ielno,i) xx(i)=x(jj) yy(i)=y(jj) do 21 j=1,8 aer(j,2*i)=0. aei(j,2*i)=0. aer(j,2*i-1)=0. aei(j,2*i-1)=0. ber(j,2*i)=0. bei(j,2*i)=0. ber(j,2*i-1)=0. bei(j,2*i-1)=0. 21 continue 20 continue ybot=y(icon(ielno,1)) yeps=0.02*(y(icon(ielno,4))-ybot) xeps=0.02*(x(icon(ielno,2))-x(icon(ielno,1))) nint=ngpint(ielno) write(*,*) ybot,yeps,xeps if(ybot.ge.1.)then c --> within buffer; very low Reynolds number, no alpha-effect remag=1. tdkill=0. else c --> within sphere remag =1. tdkill=1. endif c---------- integrate terms (gaussian quadrature of order ngp X ngp) igp=0 do 100 intx=1,nint xi=gp(nint,intx) do 90 inty=1,nint eta=gp(nint,inty) igp=igp+1 call shap23(igp) wgtjac=wt(nint,intx)*wt(nint,inty)*jac(ielno,igp) mu =xgp(ielno,igp) rr =ygp(ielno,igp) sin2=1.-mu*mu sin1=sqrt(sin2) r2 =rr*rr c --> use user-supplied functions to compute alpha, omega, and its derivatives c [tdkill is there to zero out functionals in r/R>1] alp0 = alpxy(mu,rr) *tdkill omg = omgxy(mu,rr) domgdr=(omgxy(mu,rr+yeps)-omgxy(mu,rr-yeps))/(2.*yeps) + *tdkill domgdm=(omgxy(mu+xeps,rr)-omgxy(mu-xeps,rr))/(2.*xeps) + *tdkill remag=dmagxy(mu,rr) dremag=(dmagxy(mu,rr+yeps)-dmagxy(mu,rr-yeps))/(2.*yeps) ck1 =r2 ck1ar=ck1*remag ck1am=sin2*remag ck1br=ck1*remag ck1bm=sin2*remag sigb =-1./sin2*remag c --> term multiplying dA/dmu in phi-induction eq. ckb1=-r2*sin2*domgdr *comega c --> term multiplying dA/dr in phi-induction eq. ckb3= r2*sin2*domgdm *comega c --> term multiplying A in phi-induction eq. ckb5= rr*(rr*mu*domgdr+sin2*domgdm) *comega c --> term multiplying B_y in poloidal induction equation cka1= r2*alp0 *calpha ckat= r2*tdkill do 80 i=1,4 i1=2*i-1 i2=2*i do 81 j=1,4 j1=2*j-1 j2=2*j c --> coefficients of ``poloidal'' induction equation ber(i1,j1)=ber(i1,j1)+wgtjac*ckat* phi(j) *phi(i) aer(i1,j1)=aer(i1,j1)+wgtjac*( + -ck1am *dphdx(j) *dphdx(i) + -ck1ar *dphdy(j) *dphdy(i) + +sigb *phi(j) *phi(i)) aer(i1,j2)=aer(i1,j2)+wgtjac*( + cka1 *phi(j) *phi(i)) c --> coefficients of ``toroidal'' (y) induction equation ber(i2,j2)=ber(i2,j2)+wgtjac*r2* phi(j) *phi(i) aer(i2,j1)=aer(i2,j1)+wgtjac*( + ckb1 *dphdx(j) *phi(i) + +ckb3 *dphdy(j) *phi(i) + +ckb5 *phi(j) *phi(i)) aer(i2,j2)=aer(i2,j2)+wgtjac*( + -ck1bm *dphdx(j) *dphdx(i) + -ck1br *dphdy(j) *dphdy(i) + +sigb *phi(j) *phi(i)) 81 continue 80 continue 90 continue 100 continue return c----------------------------------- end c*************************************************************************** subroutine shap23(igp) c*************************************************************************** c evaluates shape function and derivatives at points (xi,eta) * c for 4-node quadrilateral * c*************************************************************************** parameter(nm=2145,ne=2048,nb=141,ng=4) common/elvar/ xi,eta,jac(ne,ng),jacinv(ne,ng,4), + xx(4),yy(4),phi(4), + dphdx(4),dphdy(4),dphdxi(4),dphdet(4),ielno,nnd common/cntrl/ ijac real jacinv,jac c---------- constants xip=1.+xi xim=1.-xi etp=1.+eta etm=1.-eta c******************** shape functions phi(1)=xim*etm/4. phi(2)=xip*etm/4. phi(3)=xip*etp/4. phi(4)=xim*etp/4. c******************** derivatives of the shape function c---------- derivatives wrt xi at (xi,eta) dphdxi(1)=-etm/4. dphdxi(2)=etm/4. dphdxi(3)=etp/4. dphdxi(4)=-etp/4. c---------- derivatives wrt eta at (xi,eta) dphdet(1)=-xim/4. dphdet(2)=-xip/4. dphdet(3)=xip/4. dphdet(4)=xim/4. c******************** do jacobian calculations c---------- compute jacobian at (xi,eta) if(ijac.le.1)then dxdxi=0. dydet=0. dydxi=0. dxdet=0. do 210 i=1,4 dxdxi=dxdxi+dphdxi(i)*xx(i) dydet=dydet+dphdet(i)*yy(i) dydxi=dydxi+dphdxi(i)*yy(i) dxdet=dxdet+dphdet(i)*xx(i) 210 continue jac(ielno,igp)=dxdxi*dydet-dxdet*dydxi c---------- check positiveness of jacobian avg=(abs(dxdxi)+abs(dydet))/4. if(jac(ielno,igp).lt.avg*1.e-8)then write(*,*) ielno,jac(ielno,igp) stop '1 shap23' end if c---------- compute inverse of the jacobian at (xi,theta) jacinv(ielno,igp,1)=dydet/jac(ielno,igp) jacinv(ielno,igp,2)=dxdxi/jac(ielno,igp) jacinv(ielno,igp,3)=-dydxi/jac(ielno,igp) jacinv(ielno,igp,4)=-dxdet/jac(ielno,igp) end if c---------- derivatives wrt x and y at (xi,eta) do 10 i=1,4 dphdx(i)=jacinv(ielno,igp,1)*dphdxi(i) + +jacinv(ielno,igp,3)*dphdet(i) dphdy(i)=jacinv(ielno,igp,2)*dphdet(i) + +jacinv(ielno,igp,4)*dphdxi(i) 10 continue return c---------------------------------------- end c****************************************************************************** subroutine grid1(ngp,ityp) c****************************************************************************** c genere une grid d'elements quadrilateraux bilineaires * c * c --> calcul des coordonnees des noeuds et de la matrice de connectivite * c****************************************************************************** parameter(nm=2145,ne=2048,nb=141,ng=4) common/nodes/ x(nm),y(nm),numnp common/elmnts/ ntype(ne),ngpint(ne),icon(ne,4),numel common/gnode/ yy(nm-ne),xx(nm-ne),nely,nelx,nnodes nx=nelx+1 c---------- 1. calcul du "fond" do 10 i=1,nx y(i)=yy(1) x(i)=xx(i) 10 continue c---------- 2. nely couches d'elements c (nelx elements quadrilateraux bilineaires par couche) do 13 j=1,nely i2=1 do 14 i=(j-1)*nx+nx+1,j*nx+nx y(i)=yy(j+1) x(i)=xx(i2) i2=i2+1 14 continue 13 continue do 15 i=1,nely do 16 j=1,nelx iel=(i-1)*nelx+j icon(iel,1)=(i-1)*nx+j icon(iel,2)=(i-1)*nx+j+1 icon(iel,3)=i*nx+j+1 icon(iel,4)=i*nx+j ntype(iel)=ityp ngpint(iel)=ngp 16 continue 15 continue nnodes=(nely+1)*(nelx+1) numel=nelx*nely numnp=2*nnodes c write(*,1) nelx,nely,numel,nnodes,numnp 1 format(/,1x,'***** GRID:',5x,i3,' X ',i3,' Bilinear Elements',/, + 7x,'Total Elements: ',i6,/, + 7x,'Total Nodes: ',i6,/, + 7x,'Total Coefficients: ',i6) if(numel.gt.ne)then write(*,2) numel,ne stop ' *GRID*' endif if(nnodes.gt.nm)then write(*,3) numnp,nm stop ' *GRID*' endif 2 format(/,1x,'***** ERROR: numel= ',i5,5x,' ne= ',i5) 3 format(/,1x,'***** ERROR: numnp= ',i5,5x,' nm= ',i5) return end c****************************************************************************** subroutine gaussp c****************************************************************************** c calcule la position des points de gauss, ainsi que les coefficients de * c l'equation differentielle en ces points * c****************************************************************************** parameter(nm=2145,ne=2048,nb=141,ng=4) common/nodes/ x(nm),y(nm),numnp common/elmnts/ ntype(ne),ngpint(ne),icon(ne,4),numel common/gnode/ yy(nm-ne),xx(nm-ne),nely,nelx,nnodes common/phys/ xgp(ne,ng),ygp(ne,ng) sqr13=sqrt(1./3.) sqr35=sqrt(3./5.) do 1 i=1,numel ngp=ngpint(i) ymid=(y(icon(i,4))+y(icon(i,1)))/2. xmid=(x(icon(i,2))+x(icon(i,1)))/2. xl=(x(icon(i,2))-x(icon(i,1)))/2. yl=(y(icon(i,4))-y(icon(i,1)))/2. c 1 point de gauss if(ngp.eq.1)then xgp(i,1)=xmid ygp(i,1)=ymid end if c 2x2 points de gauss if(ngp.eq.2)then xgp(i,1)=xmid-xl*sqr13 xgp(i,3)=xmid+xl*sqr13 xgp(i,2)=xmid-xl*sqr13 xgp(i,4)=xmid+xl*sqr13 ygp(i,1)=ymid-yl*sqr13 ygp(i,3)=ymid-yl*sqr13 ygp(i,2)=ymid+yl*sqr13 ygp(i,4)=ymid+yl*sqr13 end if 1 continue return end c**************************************************************************** subroutine mesh2d(ipar) c**************************************************************************** c generation du maillage spatial 2-D, maillage temporel, et calcul des * c conditions limites * c**************************************************************************** parameter(nm=2145,ne=2048,nb=141,ng=4) common/nodes/ x(nm),y(nm),numnp common/elmnts/ ntype(ne),ngpint(ne),icon(ne,4),numel common/genrl/ mband,mfbw common/gnode/ yy(nm-ne),xx(nm-ne),nely,nelx,nnodes common/essbc/ uebc(2*ne),iuebc(2*ne),numebc complex uebc common/syseqs/ ka(nb,2*nm),kb(nb,2*nm),kab(nb,2*nm), + fb(2*nm),q(2*nm),bca(2*nm) complex ka,kb,kab,fb,q,bca common/cntrl/ ijac common/initq/ qsav(2*nm) complex qsav common/initc/ incr,incw,timeout,timein character*10 cdate,cclock character*64 titlea complex czero data pi/3.1415926536/ c--------------------------------------------------------------------------- nelx=32 nely=64 nel2=14 ijac=1 ityp=23 ngp=2 c---------- vertical direction (variable y, --> z ) c --> r/R=0 to 1 yi=0.0 yo=1. ny=nely-nel2+1 do 100 i=1,ny yy(i)=(yo-yi)*((float(i-1)/(ny-1))**1.0)+yi 100 continue c --> from r/R=1 to 5. yi=1. yo=5. do 101 i=1,nel2 yy(ny+i)=(yo-yi)*((float(i)/nel2)**2.0)+yi 101 continue c---------- angular direction (variable x, --> mu=cos(theta) ) c Mesh is in constant x-increments xi=-pi xo=0. nx=nelx+1 do 104 i=1,nx xx(i)=(xo-xi)*(((float(i-1)/nelx)))+xi xx(i)=cos(xx(i)) 104 continue c---------- generate grid and connectivity matrix call grid1(ngp,ityp) c---------- compute location of gauss points + some PDE coefficients call gaussp c---------- initialisations call initdt c---------- calcul de la demie-largeur de bande dU systeme if(ityp.eq.23) mband = (2*nelx+7) if(ityp.eq.24) mband = (6*nelx+11) mfbw=2*mband-1 if(mfbw.gt.nb)then write(*,5) mfbw,nb stop ' *MESH*' endif 5 format(/,1x,'***** ERROR: mfbw= ',i5,5x,' nb= ',i5) ii=40*(nelx+1) yeps=0.02*(y(ii+nelx+1)-y(ii-nelx-1)) do 119 i=ii+1,ii+nelx+1 domgdr=(omgxy(x(i),y(i)+yeps)-omgxy(x(i),y(i)-yeps))/2./yeps write(*,118) x(i),y(i),omgxy(x(i),y(i)),alpxy(x(i),y(i)),domgdr 119 continue 118 format(5f12.4) c stop ' *TMP*' c---------- guess eigenvector for very first step rzc=0.7 ybl=0.1 czero=cmplx(0.,0.) do 120 i=1,nnodes i1=2*i-1 i2=2*i theta=acos(x(i)) c --> seed toroidal field within BL if(y(i).le.1.0)then c if(y(i).le.(rzc+ybl).and.y(i).ge.(rzc-ybl))then frac=sin(pi*(y(i)-(rzc-ybl))/2./ybl ) c aq1= y(i)*(cos(theta)+sin(2.*theta)) aq1= y(i)*(1.-y(i))*(cos(theta)) aq1= y(i)*(1.-y(i))*(cos(theta)+sin(2.*theta)) c aq1= frac*(cos(theta)+sin(2.*theta)) qsav(i2)=cmplx(aq1,aq1) else qsav(i2)=czero endif c --> no xz-field (``poloidal field'') qsav(i1)=czero 120 continue c---------- conditions limites ------------------------------------------ c c 1. condition de dirichlet ("essential boundary conditions") numebc=0 c --> conditions at r=0: B_y=0, A=0 ii=numebc do 140 i=1,nelx i1=2*i-1 i2=2*i iuebc(ii+1)=i1 iuebc(ii+2)=i2 uebc (ii+1)=czero uebc (ii+2)=czero ii=ii+2 write(*,*) x(i),y(i) 140 continue numebc=ii c --> conditions within buffer zone: B_phi=0 ii=numebc do 141 i=nnodes,nnodes-((nel2+1)*(nelx+1))+1,-1 i2=2*i iuebc(ii+1)=i2 uebc (ii+1)=czero ii=ii+1 write(*,*) x(i),y(i) 141 continue numebc=ii c --> conditions along ``left'' boundary: B_y=0, A=0 ii=numebc do 142 i=1,nnodes,nelx+1 i1=2*i-1 i2=2*i iuebc(ii+1)=i1 iuebc(ii+2)=i2 uebc (ii+1)=czero uebc (ii+2)=czero ii=ii+2 write(*,*) x(i),y(i) 142 continue numebc=ii c --> conditions along ``right'' boundary: B_y=0, A=0 ii=numebc do 143 i=nelx+1,nnodes,nelx+1 i1=2*i-1 i2=2*i iuebc(ii+1)=i1 iuebc(ii+2)=i2 uebc (ii+1)=czero uebc (ii+2)=czero ii=ii+2 write(*,*) x(i),y(i) 143 continue numebc=ii c --> conditions along equatorial plane (only if ipar =1 or =2) if(ipar.gt.0)then ii=numebc do 144 i=nelx/2+1,nnodes,nelx+1 i1=2*i-1 i2=2*i if(ipar.eq.2) iuebc(ii+1)=i1 if(ipar.eq.1) iuebc(ii+1)=i2 uebc (ii+1)=czero ii=ii+1 write(*,*) x(i),y(i) 144 continue numebc=ii endif c stop ' *TMP*' return c......................................................................... end FUNCTION gammp(a,x) REAL a,gammp,x CU USES gcf,gser REAL gammcf,gamser,gln if(x.lt.0..or.a.le.0.)pause 'bad arguments in gammp' if(x.lt.a+1.)then call gser(gamser,a,x,gln) gammp=gamser else call gcf(gammcf,a,x,gln) gammp=1.-gammcf endif return END C (C) Copr. 1986-92 Numerical Recipes Software Vs1&v%1jw# following is for constant eta dmagxy=1. return end c**************************************************************************** function omgxy(xx,yy) c**************************************************************************** c computes angular velocity functional within model (0 following is for cylindrical isorotation, cf. eq. (5.10) omgxy =(yy**2) *(1.-xx**2) return end c**************************************************************************** function alpxy(xx,yy) c**************************************************************************** c computes alpha-effect functional within model (0 following is for alpha-effect as in eq. (5.11) alpxy=g1r *xx return end