[Nek5000-users] userf Questions

[nek5000-users at lists.mcs.anl.gov]

Hi All,

Hate to be so long winded, but I'm seeing problems with my code and was 
wondering if anyone could help...

I'm modeling a wind turbine using a modified actuator line method, where 
I describe the blades (3) as a set of rotating rectangular outlines, 
advanced in space at each time step.  I'm selecting GLL points which 
fall within my outlines, reading in the velocities, calculating the body 
forces, then assigning them back to ffx,ffy,ffz.  This generates a 
downstream flowfield and wake in NEK...  Looking at rough results, they 
compare (velocity and vorticity plots and isosurfaces) favorably with 
other work. However, during some recent debugging to add some features 
to my code, I noticed that when I repeatedly ran short runs and printed 
results at the same time step, my point selections for the 3 blades 
varied with each run, even though I made no other changes to the code or 
any runtime parameters.  A chunk of my point selection for a given area 
of blade passage looks like:

---------------------F77 Code In 
userf-----------------------------------------------------------------------------------
C     Test for location in Quadrant...
C
       if (thetacalc .GE. (pi/2.0 - epsim2)) then
C
C     Close to Vertical Axis...
C
C     Sin+ ...
C
       xlo = centerx + (rhubm*SIN(thetacalc))
       xhi = centerx + (bladeradius*SIN(thetacalc))

       ymax=centery+(1./SIN(thetacalc))*(chordbig/2.0)
       ymin=centery-(1./SIN(thetacalc))*(chordbig/2.0)

       if (((x.GE.xlo).AND.(x.LE.xhi)).AND.
      &     ((y.GE.ymin).AND.(y.LE.ymax))) then
C
C     Calculate forces...
C
       radm = SQRT((x-centerx)**2.0 + (y-centery)**2.0)
       vthetam = SQRT(ux**2.0 + uy**2.0)
       psim = ATAN(uz/(avelm*radm + vthetam))   ....
..........
..........
-----------------------F77 Code 
Continues------------------------------------------------------

The xlo,xhi,ymax,ymin values are calculated based on the rotational 
location of the blade, which changes at each time step, and defines an 
outline of the blade.  Then the IF statements test the x,y point (z-axis 
location has already been tested) to be within this blade outline, and 
when a point "HITS" any one of these outlines for the 3 blades, it 
calculates the forces, prints a diagnostic line, etc... A typical print 
clip, which only runs at one time step I define and prints for the "Hit" 
point, is:

---------------------F77 Code In 
userf-----------------------------------------------------------------------------------

C     Print output
C
       if(istep.eq.pstepm)then
       WRITE(4,1000)'Q1V',x,y,z,ux,uy,uz,psim,ffx,ffy,ffz
       endif
........
........
-----------------------F77 Code 
Continues------------------------------------------------------

This tells me the quadrant and general area of the point; because of 
geometric issues, I have separate selection (and associated print) 
routines for each Cartesian quadrant which are all tested in turn...

So...  My problem is if I run this code over and over, and print at say 
istep=10, the code shows a different number of points selected for each 
of the 3 blades, on each run!  Oddly, there are repeated patterns, and 
the total is always the same; for this particular setup at this istep, 
it was 204 points for all 3 blades.  An example of the differences is:

                     Blade 1                    Blade 
2                    Blade 3
-----------------------------------------------------------------------------
Run 1          76 125                           3

Run 2         32 110                           62

Run 3          78 79                            47

I've tried a bunch of things to fix it, including going back to a 
previous version which had no SAVE statements in it (had more manual 
parameter input); compiling between each run, etc...

All my code is within userf.  I'm running this in MPI mode on a large 
SMP server/workstation; compiled with gfortran.  I've run this out to 
10000 cycles with repeated runs and have always gotten qualitatively 
identical results.  I ran into this problem when I was debugging a new 
version which added a bunch of automated features I've been working on.  
But, when going back to previous versions and performing repeated short 
runs, I discovered that it did it there too...

Am I missing something here with MPI and my selection code not working 
correctly because of the way userf is called?  My understanding of the 
code is that at each time step, all the GLL points and elements are 
cycled through with a call to userf for each point; the forces are 
brought back in as ffx,ffy,ffz.  If I look qualitatively at a typical 
run through time, you can see each blade rotating step-by-step, and see 
the physical outline and edges of the blades.  It makes me think that my 
selection code (which I've debugged and checked out a lot in the past) 
is basically correct; just don't know why it would select different 
points with those IF statements on repeated runs.  Now, I can see some 
roundoff-resolution issues where it might vary one or two points between 
runs, but here the differences are huge.

Thanks for listening ;-)

Murph

Murphy Leo O'Dea
PhD Student
Mechanical Engineering
Oakland University
mlodea at oakland.edu