C----------------------------------------------------------------------- C Cicular Polarized Flow of Galloway & Proctor C (Helical Forcing Dynamo) C C Magnetic field growth rate (or 1/2 of magnetic energy growth rate) C is equal to ~ 0.3 for K_z=0.57 at Rm=100 (magnetic Reynolds #) C (Galloway & Proctor, Nature, 356, pp 692--693, 1992; C Cattaneo, Kim, Proctor & Tao, Phys. Rev. Lett., 75(8), C pp. 1523--1525, 1995 (CKPT95)) C C Note: magnetic field growth rate (or 1/2 of magnetic energy growth C rate) is equal to ~ 0.26 for K_z=0.57 at Rm=10 (CKPT95) C C Files: C gpf.* Base case (Re=Rm=10) C gpf_re.rea Case w/ Re=Rm=100 C gpf_b.* Base case for GTP solver w/ .box file C gpf_m.* Base case for GTP solver w/ .map file C C MHD: C C For an MHD case, B-field is stored after temperature array C (and passive scalars if any) so ldimt>=2. Note that in all routines C w/ NEKUSE common block i.e. uservp, userf, userq, userbc & useric, C B-field can be accessed through ux,uy,uz when C they are called with ifield=ifldmhd (ifield=1 is for velocity, C ifield=2 for temperature or ifield=ifldmhd=2 when IFHEAT=.FALSE., C etc.) in addition to the direct reference bx(ix,iy,iz,ieg), C by(ix,iy,iz,ieg), etc. C Also note that a call to outpost in userchk enforces B-field dump C in odd numbers of .f/.fld files in place of velocity. Then for C restart, one has to provide two files in .rea file -- first, C velocity .f/.fld (even number) file name and on the next line, a C file name of the B-field dump (odd .f/.fld number) both preceeded C with a line C 2 PRESOLVE/RESTART OPTIONS C C A major limitation of the current version of MHD implementation, C is the requirement that the boundary conditions for magnetic field C have to be of the same type as for velocity field, e.g., C both periodic ('PER'), Direchlet ('v '), etc. C Also ioinfo control and avg_all routine do not have a support C for magnetic field yet. C C C C GTP (Global Tensor Product Solver) : C C param(116)=p116, p117 & p118 can be used to activate global C tensor product solver based on global fast-diagonalization method C (for box-like geometries only) which alternatively, can also be C switched on by ndim<0 in .rea file on the line w/ "NEL,NDIM,NELV". C In the latter (recommended) mode of ndim<0, the mesh w/ bcs are read C from .box file, nel should be positive but no mesh w/ bcs should be C in .rea file, no .map file is needed and p116 -- p118 are ignored. C In alternative mode, one may still use .map file (and .rea/.re2 mesh w/ bcs) C with gtp/gfdm solver by setting p116, p117 & p118 to nelx, nely & nelz, C correspondingly, with an option of being multiplied by a minus sign C that indicates a primary coordinate direction for element partitioning. C If p116 -- p118 are positive or ndim<0 then the default partitioning C for gtp/gfdm solver is in nely*nelz pencils of length nelx so nely*nelz C should be greater than a number of processors np or for greater parallel C effeciency, to be a multiple of np. C Note that to use gtp/gfdm solver one also needs lelx, lely, lelz C in SIZE file to be equal or greater than nelx, nely, nelz, C correspondingly, in addition to uncommenting the second pair of C lines w/ lelg_sm & ltfdm2 in nek5_svn/trunk/nek/ZPER (make sure C you use 'make clean' first) C The solver enforces divergence-free condition upto machine accuracy C and is especilly good for meshes w/ elements of high aspect ratio C but the solver is not necessary for MHD problems with Nek5000 C so one can use standard Nek5000 solver w/ .map file by keeping C zero p116 -- p118 and positive ndim. C C----------------------------------------------------------------------- subroutine userchk include 'SIZE' include 'TOTAL' c common /chk/ energy,energy0 common /test/ uux(lx1,ly1,lz1,lelt),uuy(lx1,ly1,lz1,lelt), $ uuz(lx1,ly1,lz1,lelt) c parameter(lt=lx1*ly1*lz1*lelv) c common /test/ uux(lt),uuy(lt), c $ uuz(lt) save /test/ real unorm,unorm1,kinerg1,kinerg2 integer icalld save icalld data icalld /0/ integer d,e,m if (icalld.eq.0) then ifheat = .false. ! no heat equation to solve icalld = 1 endif m = nx1*ny1*nz1*nelv kinerg1 = 0.5*(glsc3(vx,bm1,vx,m) $ +glsc3(vy,bm1,vy,m) $ +glsc3(vz,bm1,vz,m))/volvm1 m = nx1*ny1*nz1*nelv c m = lt c if (ifield.eq.1) then ! velocity c omega = 1 c st = sin(omega*(time)) c ct = cos(omega*(time)) do i=1,m x = xm1(i,1,1,1) y = ym1(i,1,1,1) z = zm1(i,1,1,1) c R1 = sqrt(x*x+y*y); if(R1 .le. 1.0) then vx(i,1,1,1) = -y vy(i,1,1,1) = -x vz(i,1,1,1) = omega else vx(i,1,1,1) = -0.0 vy(i,1,1,1) = -0.0 vz(i,1,1,1) = 0.0 endif enddo c endif c unorm = 0.5*(glsc2(uux,uux,m) c $ +glsc2(uuy,uuy,m) c $ +glsc2(uuz,uuz,m) c $ )/volvm1 c unorm1 = 0.5*(glsc2(vx,vx,m) c $ +glsc2(vy,vy,m) c $ +glsc2(vz,vz,m) c $ )/volvm1 ibstep = iostep if (istep.gt.0.and.mod(istep,ibstep).eq.0) $ call outpost(bx,by,bz,pm,t,' ') ! odd blah.fld/.f for B-field c $ call outpost(bx,by,bz,pm,t,'mag') ! magblah.fld/.f for B-field if (ifxyo.and.istep.gt.istep) ifxyo = .false. kinerg2 = 0.5*(glsc3(vx,bm1,vx,m) $ +glsc3(vy,bm1,vy,m) $ +glsc3(vz,bm1,vz,m))/volvm1 n = nx1*ny1*nz1*nelfld(ifldmhd) ! Magnetic energy !?nelb energy0 = energy energy = glsc3(bx,bm1,bx,n) $ + glsc3(by,bm1,by,n) $ + glsc3(bz,bm1,bz,n) energy = 0.5*energy/volfld(ifldmhd) if (istep.gt.1) then ratio = energy/energy0 growth = 0.5*log(ratio)/dt ! 1/2 of magnetic energy growth rate if (nid.eq.0) write(6,1) istep,time,growth,energy, $ kinerg1,kinerg2 1 format(1i8,1p5e15.6,' growth') nstepag = param(120) ! p121 < t_avg < p120 nstepsa = param(121) if (istep.lt.nstepag) call $ runtimeavg(Emga,growth,1,nstepsa,1,'gr_Em') ! avg. growth endif return end c----------------------------------------------------------------------- subroutine uservp (ix,iy,iz,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' if (ifield.eq.1) then ! velocityc utrans = param(1) udiff = param(2) elseif (ifield.eq.ifldmhd) then ! B-field utrans = 1.0 R1 = sqrt(x*x+y*y); if(R1 .le. 2.0) then udiff= param(29) else udiff= 1000.0D0*param(29) endif endif c udiff =0. c utrans=0. return end c----------------------------------------------------------------------- subroutine userf (ix,iy,iz,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' c ffx = 0.0 ffy = 0.0 ffz = 0.0 c c if (ifield.eq.1) then ! velocity R1 = sqrt(x*x+y*y); c if(R1 .gt. 1.0) then c ffx = -ux/eta c ffy = -uy/eta c ffz = -uz/eta c else c ffx = 0.0 c ffy = 0.0 c ffz = 0.0 c endif c endif c c return end c----------------------------------------------------------------------- subroutine userq (ix,iy,iz,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' C qvol = 0.0 return end c----------------------------------------------------------------------- subroutine userbc (ix,iy,iz,iside,ieg) c NOTE ::: This subroutine MAY NOT be called by every process include 'SIZE' include 'TOTAL' include 'NEKUSE' c ux=0.0 uy=0.0 uz=0.0 c return end c----------------------------------------------------------------------- subroutine useric (ix,iy,iz,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' c if (ifield.eq.1) then ! velocity c omega = 1 R1 = sqrt(x*x+y*y); if(R1 .le. 1.0) then ux = -y uy = x uz = omega else ux = 0.0 uy = 0.0 uz = 0.0 endif c elseif (ifield.eq.ifldmhd) then ! B-field c eps = 1.e-5 eps = 1.e-7 R1 = sqrt(x*x+y*y); c if(R1 .le. 1.0) then ux = eps*sin(y)*sin(-0.39*z) ! change to ran1? uy = 0.0 uz = 0.0 c else c ux = 0.0 c uy = 0.0 c uz = 0.0 c endif endif c return end c----------------------------------------------------------------------- subroutine usrdat3 include 'SIZE' include 'TOTAL' c return end c----------------------------------------------------------------------- subroutine usrdat include 'SIZE' include 'TOTAL' c return end c----------------------------------------------------------------------- subroutine usrdat2 ! Modify geometry here include 'SIZE' include 'TOTAL' c c Rescale x,y, and z. c param(66) = 4 ! .fld -- comment out for .f param(67) = 4 ! .fld -- comment out for .f return end c----------------------------------------------------------------------- subroutine runtimeavg(ay,y,j,istep1,ipostep,s5) c c Computes, stores and for ipostep!=0 prints runtime averages c of j-quantity y (along w/ y itself unless ipostep<0) c with j + 'rtavg_' + (unique) s5 every ipostep for istep>=istep1 c include 'SIZE' ! for nid & TSTEP include 'TSTEP' ! for istep & time (or TOTAL) parameter (lymax=99) ! max of averages to store character*5 s5 ! unique string to append to 'rtavg_' character*14 s14 ! i3 + 'rtavg_' + s5 real a (lymax) ! run time average character*5 as5(lymax) ! j's strings s5 save a, as5 data a /lymax*0./, as5 /lymax*' '/ if (nid.eq.0.and.istep.ge.istep1) then ! otherwsie skip istep c--- if (istep1.lt.0) then istep1 = 0 write(6,*) 'Warning: istep1<0 in runtimeavg -- resetting to 0' endif if (as5(j).ne.s5.and.as5(j).ne.' ') then ! check for unique j & s5 write(6,10) j,s5,as5(j) call exitt endif iistep = istep - istep1 + 1 if (iistep.gt.1) then ! a_ii = (1-1/ii)*a_ii-1 + 1/ii*y_ii wy = 1./iistep wa = 1. - wy ay = wy*y + wa*a(j) ! runtime average a_ii else if (iistep.eq.1) then ! skip istep