!-----------------------------------------------------------------------

!

! user subroutines required by nek5000

!

! Parameters used by this set of subroutines:

!

!-----------------------------------------------------------------------

subroutine uservp (ix,iy,iz,ieg)

include 'SIZE'

include 'NEKUSE' ! UDIFF, UTRANS

UDIFF =0.0

UTRANS=0.0

return

end

!-----------------------------------------------------------------------

subroutine userf (ix,iy,iz,ieg)

include 'SIZE'

include 'NEKUSE' ! FF[XYZ]

include 'PARALLEL'

integer ix,iy,iz,ieg

ffx = 0.0

ffy = 0.0

ffz = 0.0

return

end

!-----------------------------------------------------------------------

subroutine userq (ix,iy,iz,ieg)

include 'SIZE'

include 'NEKUSE' ! QVOL

QVOL = 0.0

return

end

!-----------------------------------------------------------------------

subroutine userchk

implicit none

include 'SIZE' !

include 'TSTEP' ! ISTEP, lastep, time

include 'INPUT' ! IF3D, PARAM

! start framework

if (ISTEP.eq.0) call frame_start

! monitor simulation

call frame_monitor

! save/load files for full-restart

call chkpt_main

! for tripping

call stat_avg

! finalise framework

if (ISTEP.eq.NSTEPS.or.LASTEP.eq.1) then

call frame_end

endif

return

end

!-----------------------------------------------------------------------

subroutine userbc (ix,iy,iz,iside,eg)

include 'SIZE'

include 'NEKUSE' ! UX, UY, UZ, TEMP, X, Y

! velocity

ux = 0.0

uy = 0.0

uz = 0.0

return

end

!-----------------------------------------------------------------------

subroutine useric (ix,iy,iz,ieg)

include 'SIZE'

include 'NEKUSE' ! UX, UY, UZ, TEMP, Z

! argument list

integer ix,iy,iz,ieg

! local variables

real xl(LDIM), eps

real fcoeff(3) ! coefficients for random distribution

! functions

real mth_ran_dst

xl(1) = X

xl(2) = Y

xl(NDIM) = Z

eps = 1.0e-03

fcoeff(1)= 3.0e4

fcoeff(2)= -1.5e3

fcoeff(3)= 0.5e5

ux=1.0+eps*mth_ran_dst(ix,iy,iz,ieg,xl,fcoeff)

fcoeff(1)= 2.3e4

fcoeff(2)= 2.3e3

fcoeff(3)= -2.0e5

uy=eps*mth_ran_dst(ix,iy,iz,ieg,xl,fcoeff)

fcoeff(1)= 2.e4

fcoeff(2)= 1.e3

fcoeff(3)= 1.e5

uz=eps*mth_ran_dst(il,jl,kl,ieg,xl,fcoeff)

return

end

!-----------------------------------------------------------------------

subroutine usrdat

include 'SIZE'

return

end

!-----------------------------------------------------------------------

subroutine usrdat2

include 'SIZE'

include 'SOLN' ! vx,vy,vz,pr,t

integer iel

c mesh stretching

real Betax, Betay

c Hill parameters

real Lx, Ly, Lz, W, H

common /hill_param/ Lx, Ly, Lz, W, H, Betax, Betay

Betax = 2.0

Betay = 2.4

Lx = 9.

Ly = 3.035

Lz = 4.5

W = 1.929

H = 1.

c Lx

c <----------------->

c ___________________

c ^

c |

c |

c _ _ | Ly

c ^ \ / |

c H | \ / |

c v \___________/ v

c <-->

c W

c Transform box mesh to periodic hill (only for the conforming mesh!!)

call box2phill

c call gen_re2(0)

return

end

!-----------------------------------------------------------------------

subroutine usrdat3

include 'SIZE'

include 'INPUT' ! param, if3d

include 'MASS' ! volvm1

c Local variables

real Ubulk

real Betax,Betay

c Hill parameters

real Lx, Ly, Lz, W, H

common /hill_param/ Lx, Ly, Lz, W, H, Betax, Betay

c apply mass flux to drive the flow such that

c - Ubulk = velocity averaged over the whole domain

c - Ubar_inlet = velocity averaged over the inlet plane = 1.

if (if3d) then

Ubulk = (Lx*(Ly-H)*Lz)/volvm1

else

Ubulk = (Lx*(Ly-H))/volvm1

endif

if (nid.eq.0) write(6,*) 'U_bulk = ', Ubulk, ' Ubar_inlet = 1'

param(54) = -1 ! x-direction

param(55) = Ubulk ! Ubulk

return

end

c-----------------------------------------------------------------------

subroutine box2phill

implicit none

include 'SIZE'

include 'GEOM' ! {x,y,z}m1

include 'INPUT' ! param

include 'SOLN'

integer i, ntot

real Betax, Betay

real Lx, Ly, Lz, Wh, H ! x dimension, y dimension, hill half width, hill height

common /hill_param/ Lx, Ly, Lz, Wh, H, Betax ,Betay

real shift, amp

real xscale, yscale, zscale, yh, xx, yy, zz

real hill_step,hill_height,xfac,glmax,glmin

real xmin, xmax, ymin, ymax, zmin, zmax

save xmin, xmax, ymin, ymax, zmin, zmax

logical ifminmax

save ifminmax

data ifminmax /.false./

ntot = nx1*ny1*nz1*nelt

if (.not.ifminmax) then

ifminmax = .true.

xmin = glmin(xm1,ntot)

xmax = glmax(xm1,ntot)

ymin = glmin(ym1,ntot)

ymax = glmax(ym1,ntot)

if (if3d) then

zmin = glmin(zm1,ntot)

zmax = glmax(zm1,ntot)

endif

endif

C decrease resolution in the high velocity regions (increase CFL)

do i=1,ntot

xm1(i,1,1,1) = 0.5*(sinh(Betax*(xm1(i,1,1,1)-0.5))/

$ sinh(Betax*0.5) + 1.0)

enddo

c increase resolution near the wall

do i=1,ntot

ym1(i,1,1,1) = 0.5*(tanh(Betay*(2*ym1(i,1,1,1)-1.0))/

$ tanh(Betay) + 1.0)

enddo

c rescale rectangular domain [0,Lx]x[0,Ly]x[0,Lz]

xscale = Lx/(xmax-xmin)

yscale = Ly/(ymax-ymin)

do i=1,ntot

xx = xm1(i,1,1,1)

yy = ym1(i,1,1,1)

xm1(i,1,1,1) = (xx - xmin) * xscale

ym1(i,1,1,1) = (yy - ymin) * yscale

enddo

if (if3d) then

zscale = Lz/(zmax-zmin)

do i=1,ntot

zz = zm1(i,1,1,1)

zm1(i,1,1,1) = (zz - zmin) * zscale

enddo

endif

c Shift points in x

amp = 0.25

do i=1,ntot

xx = xm1(i,1,1,1)

yy = ym1(i,1,1,1)

xm1(i,1,1,1) = xx + amp*shift(xx,yy,Lx,Ly,Wh)

enddo

c Add hill

do i=1,ntot

xx = xm1(i,1,1,1)

yy = ym1(i,1,1,1)

yh = hill_height(xx,Lx,Wh,H)

yscale = 1-yh/Ly

ym1(i,1,1,1) = yh + yy * yscale

enddo

return

end

c-----------------------------------------------------------------------

C Step function for the hill

C

C x=0

C |

C _____|

C ^ \

C | \ x->

C h | \

C v \_____

C <-->

C w

function hill_step(x,w,h)

implicit none

real x,xs,w,h

real y,hill_step

xs = x/w

if (xs.le.0) then

y = h

elseif (xs.gt.0.and.xs.le.9./54.) then

y = h*min(1.,1.+7.05575248e-1*xs**2-1.1947737203e1*xs**3)

elseif (xs.gt.9./54.and.xs.le.14./54.) then

y = h*(0.895484248+1.881283544*xs-10.582126017*xs**2

$ +10.627665327*xs**3)

elseif (xs.gt.14./54.and.xs.le.20./54.) then

y = h*(0.92128609+1.582719366*xs-9.430521329*xs**2

$ +9.147030728*xs**3)

elseif (xs.gt.20./54..and.xs.le.30./54.) then

y = h*(1.445155365-2.660621763*xs+2.026499719*xs**2

$ -1.164288215*xs**3)

elseif (xs.gt.30./54..and.xs.le.40./54.) then

y = h*(0.640164762+1.6863274926*xs-5.798008941*xs**2

$ +3.530416981*xs**3)

elseif (xs.gt.40./54..and.xs.le.1.) then

y = h*(2.013932568-3.877432121*xs+1.713066537*xs**2

$ +0.150433015*xs**3)

else

y = 0.

endif

hill_step = y

return

end

c-----------------------------------------------------------------------

function hill_height(x,Lx,w,H)

implicit none

real hill_height,hill_step,x,Lx,w,h

real xx

if (x.lt.0) then

xx = Lx + mod(x, Lx)

elseif (x.gt.Lx) then

xx = mod(x, Lx)

else

xx = x

endif

hill_height = hill_step(xx,w,h) + hill_step(Lx-xx,w,h)

return

end

c-----------------------------------------------------------------------

function shift(x,y,Lx,Ly,W)

implicit none

real x,y,Lx,Ly,W

real xfac,yfac,shift

yfac = (1-y/Ly)**3

if (x.le.W/2) then

xfac = -2./W * x

elseif (x.gt.W/2.and.x.le.Lx-W/2) then

xfac = 2./(Lx-W) * x -1-W/(Lx-W)

elseif (x.gt.Lx-W/2) then

xfac = -2./W * x + 2*Lx/W

endif

shift = xfac*yfac

return

end

!======================================================================

!> @brief Register user specified modules

subroutine frame_usr_register

implicit none

include 'SIZE'

include 'FRAMELP'

!-----------------------------------------------------------------------

! register modules

call io_register

call chkpt_register

call stat_register

return

end subroutine

!======================================================================

!> @brief Initialise user specified modules

subroutine frame_usr_init

implicit none

include 'SIZE'

include 'FRAMELP'

include 'SOLN'

!-----------------------------------------------------------------------

! initialise modules

call chkpt_init

call stat_init

return

end subroutine

!======================================================================

!> @brief Finalise user specified modules

subroutine frame_usr_end

implicit none

include 'SIZE'

include 'FRAMELP'

!-----------------------------------------------------------------------

! finalise modules

return

end subroutine

!======================================================================

!> @brief Provide element coordinates and local numbers (user interface)

!! @param[out] idir mapping (uniform) direction

!! @param[out] ctrs 2D element centres

!! @param[out] cell local element numberring

!! @param[in] lctrs1,lctrs2 array sizes

!! @param[out] nelsort number of local 3D elements to sort

!! @param[out] map_xm1, map_ym1 2D coordinates of mapped elements

!! @param[out] ierr error flag

subroutine user_map2d_get(idir,ctrs,cell,lctrs1,lctrs2,nelsort,

$ map_xm1,map_ym1,ierr)

implicit none

include 'SIZE'

include 'INPUT' ! [XYZ]C

include 'GEOM' ! [XYZ]M1

! argument list

integer idir

integer lctrs1,lctrs2

real ctrs(lctrs1,lctrs2) ! 2D element centres and diagonals

integer cell(lctrs2) ! local element numberring

integer nelsort ! number of local 3D elements to sort

real map_xm1(lx1,lz1,lelt), map_ym1(lx1,lz1,lelt)

integer ierr ! error flag

! local variables

integer ntot ! tmp array size for copying

integer el ,il ,jl ! loop indexes

integer nvert ! vertex number

real rnvert ! 1/nvert

real xmid,ymid ! 2D element centre

real xmin,xmax,ymin,ymax ! to get approximate element diagonal

integer ifc ! face number

! dummy arrays

real xcoord(8,LELT), ycoord(8,LELT) ! tmp vertex coordinates

#ifdef DEBUG

! for testing

character*3 str1, str2

integer iunit, ierrl

! call number

integer icalldl

save icalldl

data icalldl /0/

#endif

!-----------------------------------------------------------------------

! initial error flag

ierr = 0

! set important parameters

! uniform direction; should be taken as input parameter

! x-> 1, y-> 2, z-> 3

idir = 3

! get element midpoints

! vertex number

nvert = 2**NDIM

rnvert= 1.0/real(nvert)

! eliminate uniform direction

ntot = 8*NELV

if (idir.EQ.1) then ! uniform X

call copy(xcoord,YC,ntot) ! copy y

call copy(ycoord,ZC,ntot) ! copy z

elseif (idir.EQ.2) then ! uniform Y

call copy(xcoord,XC,ntot) ! copy x

call copy(ycoord,ZC,ntot) ! copy z

elseif (idir.EQ.3) then ! uniform Z

call copy(xcoord,XC,ntot) ! copy x

call copy(ycoord,YC,ntot) ! copy y

endif

! set initial number of elements to sort

nelsort = 0

call izero(cell,NELT)

! for every element

do el=1,NELV

! element centre

xmid = xcoord(1,el)

ymid = ycoord(1,el)

! element diagonal

xmin = xmid

xmax = xmid

ymin = ymid

ymax = ymid

do il=2,nvert

xmid=xmid+xcoord(il,el)

ymid=ymid+ycoord(il,el)

xmin = min(xmin,xcoord(il,el))

xmax = max(xmax,xcoord(il,el))

ymin = min(ymin,ycoord(il,el))

ymax = max(ymax,ycoord(il,el))

enddo

xmid = xmid*rnvert

ymid = ymid*rnvert

! count elements to sort

nelsort = nelsort + 1

! 2D position

! in general this coud involve some curvilinear transform

ctrs(1,nelsort)=xmid

ctrs(2,nelsort)=ymid

! reference distance

ctrs(3,nelsort)=sqrt((xmax-xmin)**2 + (ymax-ymin)**2)

if (ctrs(3,nelsort).eq.0.0) then

ierr = 1

return

endif

! element index

cell(nelsort) = el

enddo

! provide 2D mesh

! in general this coud involve some curvilinear transform

if (idir.EQ.1) then ! uniform X

ifc = 4

do el=1,NELV

call ftovec(map_xm1(1,1,el),ym1,el,ifc,nx1,ny1,nz1)

call ftovec(map_ym1(1,1,el),zm1,el,ifc,nx1,ny1,nz1)

enddo

elseif (idir.eq.2) then ! uniform y

ifc = 1

do el=1,nelv

call ftovec(map_xm1(1,1,el),xm1,el,ifc,nx1,ny1,nz1)

call ftovec(map_ym1(1,1,el),zm1,el,ifc,nx1,ny1,nz1)

enddo

elseif (idir.eq.3) then ! uniform z

ifc = 5

do el=1,nelv

call ftovec(map_xm1(1,1,el),xm1,el,ifc,nx1,ny1,nz1)

call ftovec(map_ym1(1,1,el),ym1,el,ifc,nx1,ny1,nz1)

enddo

endif

#ifdef DEBUG

! testing

! to output refinement

icalldl = icalldl+1

call io_file_freeid(iunit, ierrl)

write(str1,'(i3.3)') NID

write(str2,'(i3.3)') icalldl

open(unit=iunit,file='map2d_usr.txt'//str1//'i'//str2)

write(iunit,*) idir, NELV, nelsort

write(iunit,*) 'Centre coordinates and cells'

do el=1,nelsort

write(iunit,*) el, ctrs(:,el), cell(el)

enddo

write(iunit,*) 'GLL coordinates'

do el=1,nelsort

write(iunit,*) 'Element ', el

write(iunit,*) 'XM1'

do il=1,nz1

write(iunit,*) (map_xm1(jl,il,el),jl=1,nx1)

enddo

write(iunit,*) 'YM1'

do il=1,nz1

write(iunit,*) (map_ym1(jl,il,el),jl=1,nx1)

enddo

enddo

close(iunit)

#endif

return

end subroutine

!=======================================================================

!> @brief Provide velocity, deriv. and vort. in required coordinates

!! @param[out] lvel velocity

!! @param[out] dudx,dvdx,dwdx velocity derivatives

!! @param[out] vort vorticity

subroutine user_stat_trnsv(lvel,dudx,dvdx,dwdx,vort)

implicit none

include 'SIZE'

include 'SOLN'

include 'INPUT' ! if3d

! argument list

real lvel(LX1,LY1,LZ1,LELT,3) ! velocity array

real dudx(LX1,LY1,LZ1,LELT,3) ! velocity derivatives; U

real dvdx(LX1,LY1,LZ1,LELT,3) ! V

real dwdx(LX1,LY1,LZ1,LELT,3) ! W

real vort(LX1,LY1,LZ1,LELT,3) ! vorticity

! local variables

integer itmp ! dummy variable

!-----------------------------------------------------------------------

! Velocity transformation; simple copy

itmp = NX1*NY1*NZ1*NELV

call copy(lvel(1,1,1,1,1),VX,itmp)

call copy(lvel(1,1,1,1,2),VY,itmp)

call copy(lvel(1,1,1,1,3),VZ,itmp)

! Derivative transformation

! No transformation

call gradm1(dudx(1,1,1,1,1),dudx(1,1,1,1,2),dudx(1,1,1,1,3),

$ lvel(1,1,1,1,1))

call gradm1(dvdx(1,1,1,1,1),dvdx(1,1,1,1,2),dvdx(1,1,1,1,3),

$ lvel(1,1,1,1,2))

call gradm1(dwdx(1,1,1,1,1),dwdx(1,1,1,1,2),dwdx(1,1,1,1,3),

$ lvel(1,1,1,1,3))

! get vorticity

if (IF3D) then

! curlx

call sub3(vort(1,1,1,1,1),dwdx(1,1,1,1,2),

$ dvdx(1,1,1,1,3),itmp)

! curly

call sub3(vort(1,1,1,1,2),dudx(1,1,1,1,3),

$ dwdx(1,1,1,1,1),itmp)

endif

! curlz

call sub3(vort(1,1,1,1,3),dvdx(1,1,1,1,1),dudx(1,1,1,1,2),itmp)

return

end subroutine

!======================================================================

! vim: set ft=fortran