Skip to content

GitLab

Commit dac9ddd1 authored by Cody Permann's avatar Cody Permann
Browse files

Working files

parent 24cfa7cb
Expand all Hide whitespace changes
Inline Side-by-side
examples/template_example/SIZE
View file @ dac9ddd1
......@@ -2,22 +2,22 @@ C Dimension file to be included
C
C HCUBE array dimensions
C
parameter (ldim=2)
parameter (lx1=8,ly1=lx1,lz1=1,lelt=300,lelv=lelt)
parameter (lxd=12,lyd=lxd,lzd=1)
parameter (lelx=20,lely=20,lelz=1)
parameter (ldim=3)
parameter (lx1=5,ly1=lx1,lz1=lx1,lelt=480,lelv=lelt)
parameter (lxd=8,lyd=lxd,lzd=lxd)
parameter (lelx=1,lely=1,lelz=1)
parameter (lzl=3 + 2*(ldim-3))
parameter (lx2=lx1-2)
parameter (ly2=ly1-2)
parameter (lz2=lz1 )
parameter (lz2=lz1-2)
parameter (lx3=lx1)
parameter (ly3=ly1)
parameter (lz3=lz1)
parameter (lp = 512)
parameter (lelg = 4100)
parameter (lp = 32)
parameter (lelg = 500)
c
c parameter (lpelv=lelv,lpelt=lelt,lpert=3) ! perturbation
c parameter (lpx1=lx1,lpy1=ly1,lpz1=lz1) ! array sizes
......@@ -36,8 +36,8 @@ c
parameter (lbx2=1,lby2=1,lbz2=1)
C LX1M=LX1 when there are moving meshes; =1 otherwise
parameter (lx1m=1,ly1m=1,lz1m=1)
parameter (ldimt= 2) ! 2 passive scalars + T
parameter (lx1m=lx1,ly1m=ly1,lz1m=lz1)
parameter (ldimt= 4) ! 3 passive scalars + T
parameter (ldimt1=ldimt+1)
parameter (ldimt3=ldimt+3)
c
......@@ -69,8 +69,8 @@ C
C
C Uzawa projection array dimensions
C
parameter (mxprev = 20)
parameter (lgmres = 30)
parameter (mxprev = 40)
parameter (lgmres = 40)
C
C Split projection array dimensions
C
......@@ -92,11 +92,11 @@ c automatically added by makenek
parameter(lxo = lx1) ! max output grid size (lxo>=lx1)
c automatically added by makenek
parameter(lpart = 1 ) ! max number of particles
parameter(lpart = 62000) ! max number of particles
c automatically added by makenek
integer ax1,ay1,az1,ax2,ay2,az2
parameter (ax1=lx1,ay1=ly1,az1=lz1,ax2=lx2,ay2=ly2,az2=lz2) ! running averages
parameter (ax1=1,ay1=1,az1=1,ax2=1,ay2=1,az2=1) ! running averages
c automatically added by makenek
parameter (lxs=1,lys=lxs,lzs=(lxs-1)*(ldim-2)+1) !New Pressure Preconditioner
......
examples/template_example/template_example.usr
View file @ dac9ddd1
c-----------------------------------------------------------------------
subroutine exact(uu,vv,xx,yy,n,time,visc,u0,v0)
c
c This routine creates initial conditions for an exact solution
c to the Navier-Stokes equations based on the paper of Walsh [1],
c with an additional translational velocity (u0,v0).
c
c The computational domain is [0,2pi]^2 with doubly-periodic
c boundary conditions.
c
c Walsh's solution consists of an array of vortices determined
c as a linear combinations of eigenfunctions of having form:
c
c cos(pi m x)cos(pi n y), cos(pi m x)sin(pi n y)
c sin(pi m x)cos(pi n y), sin(pi m x)sin(pi n y)
c
c and
c
c cos(pi k x)cos(pi l y), cos(pi k x)sin(pi l y)
c sin(pi k x)cos(pi l y), sin(pi k x)sin(pi l y)
c
c While there are constraints on admissible (m,n),(k,l)
c pairings, Walsh shows that there is a large class of
c possible pairings that give rise to very complex vortex
c patterns.
c
c Walsh's solution applies either to unsteady Stokes or
c unsteady Navier-Stokes. The solution is a non-translating
c decaying array of vortices that decays at the rate
c
c exp ( -4 pi^2 (m^2+n^2) visc time ),
c
c with (m^2+n^2) = (k^2+l^2). A nearly stationary state may
c be obtained by taking the viscosity to be extremely small,
c so the effective decay is negligible. This limit, however,
c leads to an unstable state, thus diminsishing the value of
c Walsh's solution as a high-Reynolds number test case.
c
c It is possible to extend Walsh's solution to a stable convectively-
c dominated case by simulating an array of vortices that translate
c at arbitrary speed by adding a constant to the initial velocity field.
c This approach provides a good test for convection-diffusion dynamics.
c
c The approach can also be extended to incompressible MHD with unit
c magnetic Prandtl number Pm.
c
c [1] Owen Walsh, "Eddy Solutions of the Navier-Stokes Equations,"
c in The Navier-Stokes Equations II - Theory and Numerical Methods,
c Proceedings, Oberwolfach 1991, J.G. Heywood, K. Masuda,
c R. Rautmann, S.A. Solonnikov, Eds., Springer-Verlag, pp. 306--309
c (1992).
c
c 2/23/02; 6/2/09; pff
c
c
include 'SIZE'
include 'INPUT'
c
real uu(n),vv(n),xx(n),yy(n)
c
real cpsi(2,5), a(2,5)
save cpsi , a
c data a / .4,.45 , .4,.2 , -.2,-.1 , .2,.05, -.09,-.1 / ! See eddy.m
c data cpsi / 0,65 , 16,63 , 25,60 , 33,56 , 39,52 / ! See squares.f
c data cpsi / 0,85 , 13,84 , 36,77 , 40,75 , 51,68 /
c This data from Walsh's Figure 1 [1]:
data a / -.2,-.2, .25,0., 0,0 , 0,0 , 0,0 /
data cpsi / 0, 5 , 3, 4 , 0,0 , 0,0 , 0,0 /
one = 1.
pi = 4.*atan(one)
aa = cpsi(2,1)**2
arg = -visc*time*aa ! domain is [0:2pi]
e = exp(arg)
c
c ux = psi_y, uy = -psi_x
c
do i=1,n
x = xx(i) - u0*time
y = yy(i) - v0*time
sx = sin(cpsi(2,1)*x)
cx = cos(cpsi(2,1)*x)
sy = sin(cpsi(2,1)*y)
cy = cos(cpsi(2,1)*y)
u = a(1,1)*cpsi(2,1)*cy
v = a(2,1)*cpsi(2,1)*sx
do k=2,5
s1x = sin(cpsi(1,k)*x)
c1x = cos(cpsi(1,k)*x)
s2x = sin(cpsi(2,k)*x)
c2x = cos(cpsi(2,k)*x)
s1y = sin(cpsi(1,k)*y)
c1y = cos(cpsi(1,k)*y)
s2y = sin(cpsi(2,k)*y)
c2y = cos(cpsi(2,k)*y)
c1 = cpsi(1,k)
c2 = cpsi(2,k)
if (k.eq.2) u = u + a(1,k)*s1x*c2y*c2
if (k.eq.2) v = v - a(1,k)*c1x*s2y*c1
if (k.eq.2) u = u - a(2,k)*s2x*c1y*c1
if (k.eq.2) v = v + a(2,k)*c2x*s1y*c2
if (k.eq.3) u = u - a(1,k)*s1x*c2y*c2
if (k.eq.3) v = v + a(1,k)*c1x*s2y*c1
if (k.eq.3) u = u - a(2,k)*c2x*c1y*c1
if (k.eq.3) v = v - a(2,k)*s2x*s1y*c2
if (k.eq.4) u = u + a(1,k)*c1x*c2y*c2
if (k.eq.4) v = v + a(1,k)*s1x*s2y*c1
if (k.eq.4) u = u + a(2,k)*c2x*c1y*c1
if (k.eq.4) v = v + a(2,k)*s2x*s1y*c2
if (k.eq.5) u = u - a(1,k)*s1x*c2y*c2
if (k.eq.5) v = v + a(1,k)*c1x*s2y*c1
if (k.eq.5) u = u - a(2,k)*s2x*c1y*c1
if (k.eq.5) v = v + a(2,k)*c2x*s1y*c2
enddo
uu(i) = u*e + u0
vv(i) = v*e + v0
enddo
return
end
C
C USER SPECIFIED ROUTINES:
C
C - boundary conditions
C - initial conditions
C - variable properties
C - local acceleration for fluid (a)
C - forcing function for passive scalar (q)
C - general purpose routine for checking errors etc.
C
c-----------------------------------------------------------------------
subroutine uservp (ix,iy,iz,ieg)
subroutine uservp (ix,iy,iz,eg)
include 'SIZE'
include 'TOTAL'
include 'NEKUSE'
C
integer e,f,eg
udiff =0.
utrans=0.
return
end
c-----------------------------------------------------------------------
subroutine userf (ix,iy,iz,ieg)
subroutine userf (ix,iy,iz,eg)
include 'SIZE'
include 'TOTAL'
include 'NEKUSE'
C
ffx = 0.0
ffy = 0.0
ffz = 0.0
integer e,f,eg
ffx=0.0
ffy=0.0
ffz=0.0
return
end
c-----------------------------------------------------------------------
subroutine userq (ix,iy,iz,ieg)
subroutine userq (ix,iy,iz,eg)
include 'SIZE'
include 'TOTAL'
include 'NEKUSE'
C
integer e,f,eg
c e = gllel(eg)
qvol = 0.0
source = 0.0
return
end
c-----------------------------------------------------------------------
subroutine userchk
include 'SIZE'
include 'TOTAL'
c
common /exacu/ ue(lx1,ly1,lz1,lelt),ve(lx1,ly1,lz1,lelt)
common /exacd/ ud(lx1,ly1,lz1,lelt),vd(lx1,ly1,lz1,lelt)
include 'SIZE'
include 'TOTAL'
common/test_passing/ flux_moose, temp_nek
real sint, sint1, sarea, sarea1, wtmp
integer e, f
ifield = 1 ! for outpost
n1=nelt*lx1*ly1*lz1
n2=nelt*lx2*ly2*lz2
n = nx1*ny1*nz1*nelv
visc = param(2)
u0 = param(96)
v0 = param(97)
call exact (ue,ve,xm1,ym1,n,time,visc,u0,v0)
if (istep.eq.0 ) call outpost(ue,ve,vx,pr,t,' ')
pmin=glmin(pr,n2)
pmax=glmax(pr,n2)
call sub3 (ud,ue,vx,n)
call sub3 (vd,ve,vy,n)
if (istep.eq.nsteps) call outpost(ud,vd,vx,pr,t,'err')
wmin=glmin(vz,n1)
wmax=glmax(vz,n1)
umx = glamax(vx,n)
vmx = glamax(vy,n)
uex = glamax(ue,n)
vex = glamax(ve,n)
udx = glamax(ud,n)
vdx = glamax(vd,n)
tmin=glmin(t,n1)
tmax=glmax(t,n1)
if (nid.eq.0) then
write(6,11) istep,time,udx,umx,uex,u0,' X err'
write(6,11) istep,time,vdx,vmx,vex,v0,' Y err'
11 format(i5,1p5e14.6,a7)
endif
ifflow=.false.
sint1=0.0
sarea1=0.0
do e=1,lelt
do f=1,6
call surface_int(sint,sarea,t,e,f)
if (cbc(f,e,1).eq.'W ') then
sint1=sint1+sint
sarea1=sarea1+sarea
endif
enddo
enddo
call gop(sint1,wtmp,'+ ',1)
call gop(sarea1,wtmp,'+ ',1)
if (istep.le.5) then ! Reset velocity to eliminate
call copy (vx,ue,n) ! start-up contributions to
call copy (vy,ve,n) ! temporal-accuracy behavior.
temp_nek=sint1/sarea1
if (nid.eq.0) then
write(6,*)"*** Temperature: ",tmin," - ",tmax
write(6,*)"*** Av. Temperature: ",temp_nek
endif
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'
ux=0.0
uy=0.0
uz=0.0
temp=0.0
common/test_passing/ flux_moose, temp_nek
integer e,ieg
real ucx, ucy, ucz, ucy_e, yy
e=gllel(ieg)
zz = zm1(ix,iy,iz,e)
xx = xm1(ix,iy,iz,e)
yy = ym1(ix,iy,iz,e)
rr = sqrt(xx**2 + yy**2)
ux = 0.0
uy = 0.0
uz = 1.0
temp = 0.0
flux = 1.0 !flux_moose
return
end
c-----------------------------------------------------------------------
......@@ -223,26 +123,16 @@ c-----------------------------------------------------------------------
include 'TOTAL'
include 'NEKUSE'
common /exacu/ ue(lx1,ly1,lz1,lelt),ve(lx1,ly1,lz1,lelt)
common /exacd/ ud(lx1,ly1,lz1,lelt),vd(lx1,ly1,lz1,lelt)
integer icalld
save icalld
data icalld /0/
integer idum, e
save idum
data idum / 0 /
n = nx1*ny1*nz1*nelv
if (icalld.eq.0) then
icalld = icalld + 1
time = 0.
u0 = param(96)
v0 = param(97)
call exact (ue,ve,xm1,ym1,n,time,visc,u0,v0)
endif
if (idum.eq.0) idum = 99 + nid
eps = .35
ie = gllel(ieg)
ux=ue(ix,iy,iz,ie)
uy=ve(ix,iy,iz,ie)
uz=0.0
uy=0.0
ux=0.0
temp=0
return
......@@ -251,29 +141,109 @@ c-----------------------------------------------------------------------
subroutine usrdat
include 'SIZE'
include 'TOTAL'
integer e
one = 1.
twopi = 8.*atan(one)
do e=1,nelv ! Rescale mesh to [0,2pi]^2
do i=1,4 ! Assumes original domain in .rea file on [0,1]
xc(i,e) = twopi*xc(i,e)
yc(i,e) = twopi*yc(i,e)
enddo
enddo
c param(66) = 4
c param(67) = 4
return
end
c-----------------------------------------------------------------------
subroutine usrdat2
include 'SIZE'
include 'TOTAL'
param(66) = 4. ! These give the std nek binary i/o and are
param(67) = 4. ! good default values
ifuservp=.false.
return
end
c-----------------------------------------------------------------------
subroutine usrdat3
include 'SIZE'
include 'TOTAL'
c
return
end
C=======================================================================
subroutine nek_init_step()
include 'SIZE'
include 'TSTEP'
include 'INPUT'
include 'CTIMER'
real*4 papi_mflops
integer*8 papi_flops
integer icall, kstep, i, pstep
common /cht_coupler/ pstep
save icall
if (icall.eq.0) then
call nekgsync()
if (instep.eq.0) then
if(nid.eq.0) write(6,'(/,A,/,A,/)')
& ' nsteps=0 -> skip time loop',
& ' running solver in post processing mode'
else
if(nio.eq.0) write(6,'(/,A,/)') 'Starting time loop ...'
endif
isyc = 0
itime = 0
if(ifsync) isyc=1
itime = 1
call nek_comm_settings(isyc,itime)
call nek_comm_startstat()
istep = 0
endif
istep=istep+1
if (lastep .eq. 1) then
pstep=2
else
call nek_advance
pstep=2
endif
return
end
C=======================================================================
subroutine nek_step()
include 'SIZE'
include 'TSTEP'
include 'INPUT'
include 'CTIMER'
common /cht_coupler/ pstep
integer pstep
pstep=pstep+1
call heat(pstep)
return
end
C=======================================================================
subroutine nek_finalize_step()
include 'SIZE'
include 'TSTEP'
include 'INPUT'
include 'CTIMER'
common /cht_coupler/ pstep
integer pstep
real*4 papi_mflops
integer*8 papi_flops
integer icall, kstep, knstep, i
if (param(103).gt.0) call q_filter (param(103))
call setup_convect (pstep) ! Save convective velocity _after_ filter
call userchk
call prepost (.false.,'his')
call in_situ_check()
if (mod(istep,nsteps).eq.0) lastep=1
call nek_comm_settings(isyc,0)
call comment
return
end
c-----------------------------------------------------------------------
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment