Adding ethier case

ethier/ethier.par

@@ -0,0 +1,41 @@
+[GENERAL]
+#verbose = true 
+polynomialOrder = 9
+#startFrom = "restart.fld"
+stopAt = numSteps
+numSteps = 100 
+dt = 2e-03
+timeStepper = tombo3
+writeControl = simulationTime 
+writeInterval = 0.1
+
+[PRESSURE]
+residualTol = 1e-08
+
+[VELOCITY]
+boundaryTypeMap = inlet 
+residualTol = 1e-12
+density = 1.0
+viscosity = -100
+
+[TEMPERATURE]
+boundaryTypeMap = inlet 
+residualTol = 1e-12
+rhoCp = 1.0
+conductivity = -100
+
+[SCALAR01]
+boundaryTypeMap = flux
+residualTol = 1e-12
+rho = 1.0
+diffusivity = -100
+
+[CASEDATA]
+P_U0 = 0.5
+P_V0 = 0.1
+P_W0 = 0.2
+P_A0 = 0.025
+P_D0 = 0.5
+P_OMEGA = 15.0
+P_AMP = 1.5
+

ethier/ethier.usr

@@ -0,0 +1,252 @@
+c
+c  Exact fully 3D Navier-Stokes benchmark
+c   
+c  Reference:
+c
+c    C Ross Ethier, David Steinman,
+c    Exact fully 3D Navier-Stokes solutions for benchmarking,
+c    International Journal for Numerical Methods in Fluids,
+c    Volume 19, Number 5, March 1994, pages 369-375.
+c
+c  Setup:
+c
+c    [-1,1] cube centered at (0,0,0), 
+c    from t = 0 to t = 0.1, with parameters a = PI/4 and d = PI/2,
+c    and with Dirichlet boundary conditions on all faces of the cube.
+
+C-----------------------------------------------------------------------
+      subroutine userchk()
+      include 'SIZE'
+      include 'TOTAL'
+
+      real uxe  (lx1*ly1*lz1*lelv)
+      real uye  (lx1*ly1*lz1*lelv)
+      real uze  (lx1*ly1*lz1*lelv)
+      real uxerr(lx1*ly1*lz1*lelv)
+      real uyerr(lx1*ly1*lz1*lelv)
+      real uzerr(lx1*ly1*lz1*lelv)
+      real terr (lx1*ly1*lz1*lelv)
+
+      real pre  (lx2,ly2,lz2,lelv)
+      real prerr(lx2,ly2,lz2,lelv)
+
+      real err(4)
+      save err 
+
+      common /SCRNS/ wo1(lx1,ly1,lz1,lelv)
+     &              ,wo2(lx1,ly1,lz1,lelv)
+     &              ,omg(lx1*ly1*lz1*lelv,ldim)
+
+
+      COMMON /NRSSCPTR/ nrs_scptr(1)
+      integer*8         nrs_scptr
+
+      nrs_scptr(1) = loc(err)
+
+      n  = nelv*nx1*ny1*nz1
+      n2 = nelv*nx2*ny2*nz2
+
+      call exactu(uxe,uye,uze,xm1,ym1,zm1,n)
+      call exactp(pre,xm2,ym2,zm2,n2)
+
+      pbar  = glsc2(pr ,bm2,n2)/volvm2
+      pbre  = glsc2(pre,bm2,n2)/volvm2
+      pbarm = -pbre
+      call cadd(pre,pbarm,n2)     !  Make sure pr and pre have same
+      call cadd(pre,pbar,n2)      !  average pressure before comparing.
+
+      call sub3(uxerr,vx,uxe,n)
+      call sub3(uyerr,vy,uye,n)
+      call sub3(uzerr,vz,uze,n)
+      uxerrl2 = glsc3(uxerr,bm1,uxerr,n)
+      uxerrl2 = sqrt(uxerrl2)
+
+      call sub3(terr,t,uxe,n)
+      terrl2 = glsc3(terr,bm1,terr,n)
+      terrl2 = sqrt(terrl2)
+
+      call sub3(terr,t(1,1,1,1,2),uxe,n)
+      serrl2 = glsc3(terr,bm1,terr,n)
+      serrl2 = sqrt(serrl2)
+
+      call sub3(prerr,pr,pre,n2)
+      prerrl2 = glsc3(prerr,bm2,prerr,n2)
+      prerrl2 = sqrt(prerrl2)
+
+      err(1) = uxerrl2
+      err(2) = prerrl2
+      err(3) = terrl2
+      err(4) = serrl2
+
+      if (nid.eq.nio) write(6,*) istep,time,uxerrl2,prerrl2,
+     &                           terrl2,serrl2,' L2 err'
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine usrdat()   ! This routine to modify element vertices
+      include 'SIZE'
+      include 'TOTAL'
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine usrdat2()  ! This routine to modify mesh coordinates
+      include 'SIZE'
+      include 'TOTAL'
+      include 'CASEDATA'
+
+      common /scnrs/ sc_nrs(10)
+      real               sc_nrs
+
+      P_U0 = sc_nrs(1)
+      P_V0 = sc_nrs(2)
+      P_W0 = sc_nrs(3)
+
+      P_A0 = sc_nrs(4)
+      P_D0 = sc_nrs(5)
+
+      P_OMEGA = sc_nrs(6)
+      P_AMP = sc_nrs(7)
+
+      call rescale_x(xm1,-1.0,1.0)
+      call rescale_x(ym1,-1.0,1.0)
+      call rescale_x(zm1,-1.0,1.0)
+
+      do iel=1,nelt
+      do ifc=1,2*ndim
+         if (cbc(ifc,iel,1) .eq. 'EXO') boundaryID(ifc,iel) = 1
+      enddo
+      enddo
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine usrdat3()
+      include 'SIZE'
+      include 'TOTAL'
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine useric(ix,iy,iz,eg) ! set up initial conditions
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKUSE'
+      integer e,eg
+
+      call exactu(ux,uy,uz,x,y,z,1)
+      temp = ux
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine exactu(uxe,uye,uze,x,y,z,n)
+      include 'SIZE'
+      include 'TOTAL'
+      include 'CASEDATA'
+
+      real uxe(n),uye(n),uze(n)
+      real x(n),y(n),z(n)
+
+
+      real w1(lx1*ly1*lz1*lelt)
+
+      a = P_A0*pi
+      d = P_D0*pi
+
+      call uvwp_ethier(param(2),a,d,n,x,y,z,time,uxe,uye,uze,w1)
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine exactp(pre,x,y,z,n)
+      include 'SIZE'
+      include 'TOTAL'
+      include 'CASEDATA'
+
+      real pre(n)
+      real x(n),y(n),z(n)
+
+      real w1(lx1*ly1*lz1*lelt)
+      a = P_A0*pi
+      d = P_D0*pi
+
+      call uvwp_ethier(param(2),a,d,n,x,y,z,time,w1,w1,w1,pre)
+
+      call ortho(pre)
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine uvwp_ethier (nu,a,d,n,x,y,z,t,u,v,w,p)
+
+      implicit none
+      include 'CASEDATA'
+
+      integer n
+
+      real nu
+      real a
+      real cxy(n)
+      real cyz(n)
+      real czx(n)
+      real d
+      real e2t
+      real ex(n)
+      real exy(n)
+      real ey(n)
+      real eyz(n)
+      real ez(n)
+      real ezx(n)
+      integer i
+      real p(n)
+      real sxy(n)
+      real syz(n)
+      real szx(n)
+      real t
+      real u(n)
+      real v(n)
+      real w(n)
+      real x(n)
+      real y(n)
+      real z(n)
+
+      real xx, yy, zz   
+      do i = 1, n
+
+        xx = x(i) - P_U0*t
+        yy = y(i) - P_V0*t
+        zz = z(i) - P_W0*t
+
+        ex(i) = exp ( a * xx )
+        ey(i) = exp ( a * yy )
+        ez(i) = exp ( a * zz )
+
+        e2t = exp  ( - nu * d * d * t )
+
+        exy(i) = exp ( a * ( xx + yy ) )
+        eyz(i) = exp ( a * ( yy + zz ) )
+        ezx(i) = exp ( a * ( zz + xx ) )
+
+        sxy(i) = sin ( a * xx + d * yy )
+        syz(i) = sin ( a * yy + d * zz )
+        szx(i) = sin ( a * zz + d * xx )
+
+        cxy(i) = cos ( a * xx + d * yy )
+        cyz(i) = cos ( a * yy + d * zz )
+        czx(i) = cos ( a * zz + d * xx )
+
+        u(i) = - a * ( ex(i) * syz(i) + ez(i) * cxy(i) ) * e2t + P_U0
+        v(i) = - a * ( ey(i) * szx(i) + ex(i) * cyz(i) ) * e2t + P_V0
+        w(i) = - a * ( ez(i) * sxy(i) + ey(i) * czx(i) ) * e2t + P_W0
+        p(i) = - 0.5D+00 * a * a * e2t * e2t * (
+     &    + ex(i) * ex(i) + 2.0D+00 * sxy(i) * czx(i) * eyz(i)
+     &    + ey(i) * ey(i) + 2.0D+00 * syz(i) * cxy(i) * ezx(i)
+     &    + ez(i) * ez(i) + 2.0D+00 * szx(i) * cyz(i) * exy(i) )
+
+      end do
+
+      return
+      end
+c-----------------------------------------------------------------------

ethier/mode2/CASEDATA

@@ -0,0 +1,7 @@
+      real p_u0, p_v0, p_w0,
+     $ p_a0, p_d0,
+     $ p_omega, p_amp
+
+      common /casevars/ p_u0, p_v0, p_w0,
+     $ p_a0, p_d0,
+     $ p_omega, p_amp