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 real kinerg1,kinerg2 integer icalld save icalld data icalld /0/ integer d,e m = nx1*ny1*nz1*nelv kinerg1 = 0.5*(glsc3(vx,bm1,vx,m) $ +glsc3(vy,bm1,vy,m) $ +glsc3(vz,bm1,vz,m))/volvm1 if (icalld.eq.0) then ifheat = .false. ! no heat equation to solve icalld = 1 endif if (ifield.eq.1) then ! velocity omega = 1. do i=1,nx1*ny1*nz1*nelv x = xm1(i,1,1,1); y = ym1(i,1,1,1); z = zm1(i,1,1,1) st = sin(omega*time) ct = cos(omega*time) c argxc = x + ct argyc = y + ct argxs = x + st argys = y + st c visc = param(2) s = 3./2. s = sqrt(s) c cxc = s*cos(argxc) cys = s*cos(argys) sxc = s*sin(argxc) sys = s*sin(argys) c vx(i,1,1,1) = -sys vy(i,1,1,1) = -cxc vz(i,1,1,1) = sxc + cys enddo endif 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(i8,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' C udiff =0. 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 if (ifield.eq.1) then ! velocity omega = 1. st = sin(omega*time) ct = cos(omega*time) c argxc = x + ct argyc = y + ct argxs = x + st argys = y + st c visc = param(2) s = 3./2. s = sqrt(s) c cxc = s*cos(argxc) cyc = s*cos(argyc) cxs = s*cos(argxs) cys = s*cos(argys) sxc = s*sin(argxc) syc = s*sin(argyc) sxs = s*sin(argxs) sys = s*sin(argys) ffx = 0.0 ffy = 0.0 ffz = 0.0 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 st = sin(omega*time) ct = cos(omega*time) c argxc = x + ct argyc = y + ct argxs = x + st argys = y + st c visc = param(2) s = 3./2. s = sqrt(s) c cxc = s*cos(argxc) cys = s*cos(argys) sxc = s*sin(argxc) sys = s*sin(argys) c ux = -sys uy = -cxc uz = sxc + cys c elseif (ifield.eq.ifldmhd) then ! B-field c eps = 1.e-5 eps = 1.e-6 ux = eps*sin(y)*sin(0.57*z) ! change to ran1? uy = 0. uz = 0. 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 one = 1. pi = 4.*atan(one) twopi = 2.*pi zleng = twopi / 0.57 ! dominant wavelenth & wavenumber c call rescale_x(xm1,0.,twopi) call rescale_x(ym1,0.,twopi) call rescale_x(zm1,0.,zleng) 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