[Nek5000-users] Body Force And USR Questions

Wed Aug 28 20:43:54 CDT 2013

Murph,

1) Elemental arrays (such as velocity, spatial coordinates, etc) are indeed numbered by elements.  These arrays have the form (using the x velocity, for example)

vx( lx1, ly1, lz1, lelv),

where lx1, ly1, and lz1 are the number of GLL points in each master coordinate direction (r, s, and t, respectively), and lelv is the *allocated* number of local elements on each processor (nelv is the actual number of local elements).

For example, if you wanted the x velocity at the (5,4,6) nodal point for local element 10, you would access it as

vx(5,4,6,10).

2) The physical coordinates on mesh 1 (velocity grid) are stored globally in the arrays xm1, ym1 and zm1, with the same structure as "vx" above.  The x, y, and z variables you might see in .usr file examples are temporary variables that are given to the user for convenience.  For example, the subroutine userf has the form

subroutine userf(ix,iy,iz,ieg)

where ix, iy, iz are the r,s,t indices and ieg is the global element number. The call to userf would appear in a loop as something like

do iel=1,nelv

  ieg = lglel(iel)

  do iz=1,nz1
  do iy=1,ny1
  do ix=1,nx1

    x = xm1(ix,iy,iz,iel)
   ….
   call userf(ix,iy,iz,ieg)

  … (use the user-supplied values for ffx, ffy, and ffz)

  enddo
  enddo
  enddo

enddo

Note that ieg is the "global" element number (out of the total number of elements in your problem), while iel is the "local" element number (out of the number of elements on a given processor).  The process is similar for userbc and useric.

There are only physical coordinates in Nek.  The arrays xc, yc and zc store the element corner coordinates (4 points in 2D, 8 points in 3D, for each element); basically, these arrays store what you would see in the .rea file. 

3) ux, uy, and uz are temporary variables in the same vein as x, y, and z, as listed above.  The velocity arrays are vx, vy, and vz for the x, y, and z velocities, and they have the array structure as the spatial coordinates.

4) The subroutines in the .usr file are indeed only at the GLL points.  If you wanted to apply a forcing function at a specific point/line (as I'm assuming you want to do when you specify z == 71), you would need to project that to the GLL points.  Unfortunately, a projection such as this is not trivial, and I'm not aware of any routines that are already built to do something like that.  My suggestion  would be to avoid something like z==71, and instead do something like 

  if ( z.gt.(71. - eps).and.z.lt.(71. + eps)) …

where "eps" is chosen larger enough so that you know at least *some* GLL points are encapsulated in this range (so, you might need to play with this number a bit).  Also, you may want to avoid using logical symbols such as >= and == and instead use the Fortran 77 equivalents, such as ".ge." and ".eq." .  Since the base NEK code is written in F77, using F90-type relation operators *potentially* could cause some problems. 

5) userbc and userf are indeed called every time step.  "istep" is the time step, and "time" is the simulation time, not wall-clock time.  I'm not sure what itime represents.

I hope this helps get you started!

Josh
Ph.D., Texas A&M University '13

On Aug 27, 2013, at 8:11 PM, nek5000-users at lists.mcs.anl.gov wrote:

> Hi All,
> 
> I hate to ask such lengthy and simple questions, but I'm getting stuck and could use a sharp kick...  I read through source code, examples, and the docs, but could still use some help...  I have some background with SEM study and my own trivial code, so I have some understanding of the spatial discretization part of it.  I have a trivial "template" problem setup and running.
> 
> My problem is understanding what some of the variables are (in context of their use in the code), and how to use them in the .usr file to do what I want.  I'm trying to simulate a horizontal wind turbine wake flow, using the actuator disc method and then the actuator line method.  I can get some flow results by specifying constant body forces in .usr.  I have also tried a simple use of what I think is the z-velocity in an equation to calculate a body force, and this seems to work also...
> 
> I'd like to move on to more complex body force calculations, but need to know more about what some of the variables are; how they're calculated, and how they are used/called in .usr.  The following is my first pass at a list of questions:
> 
> 1) Spatially, which variables specify the actual GLL point numbers (not physical coordinates, but point count in each coordinate direction)?  Are these ix,iy,iz?  I assume they are numbered per element?
> 
> 2) Are x,y,z the real physical coordinates in space?  Are they the actual domain size?  Are there non-dimensional coordinates also?  What are the xc,yc,zc variables?
> 
> 3) Are ux, uy, uz the fluid velocities in the coordinate directions?
> 
> 4) Can the velocities and forces be specified or referred to at any physical location, or only at the GLL points?  In other words, is there internal extrapolation going on, or do I need to specify and read my forces at GLL points only?  All the subroutines seem to call based on ix,iy,iz; when I run the following, it does calculate some forces:
> 
> subroutine userf  (ix,iy,iz,ieg)
>       include 'SIZE'
>       include 'TOTAL'
>       include 'NEKUSE'
> c
>       ffx = 0.
>       ffy = 0.
>       ffz = 0.
>       
>       if ((z>70).AND.(z<71.).AND.(x>60.).AND.(x<90.)
>      +.AND.(y>60.).AND.(y<90.)) then
>       
>       ffz = -(.5*1.2*(uz**2.)*(4.*.9*(1.-.9)))
> c      ffz = (-1000.)
>       endif
> 
> c      write (6,*) ffz,ieg
> 
>       return
>       end
> 
> However, if I change my spatial selection to:
> 
> if ((z==71.).AND.(x>60.).AND.(x<90.)
>      +.AND.(y>60.).AND.(y<90.)) then
> 
> I get no force output.  If everything is at the GLL locations, are there available subroutines to extrapolate to/from physical locations? 
> 
> 5) I assume that the subroutines in .usr are called at each time step?  Is "istep" the number of the time step?  Is "time" the actual time integration time, or the wall clock?  What is "itime"?  Just want to be able to calculate physical location of a line moving on the grid based on time step...
> 
> Thanks,
> Murph
> 
> 
> -- 
> Murphy Leo O’Dea, PE    “Murph”  
> PhD Student, Mechanical Engineering
> 
> School Of Engineering And Computer Science
> Oakland University
> 
> mlodea at oakland.edu
> 
> http://www.secs.oakland.edu/~mlodea/
> 
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20130828/6984832c/attachment-0001.html>