TEST CASE results are below. I used a contour interval of 1.0 for all of the contour plots.
Click on any image to jump to the web page for that scheme; click on ">" to step forward.

Lax-Wendroff	6th-order Crowley	Takacs

FAQs, Questions & Answers

1. [ Problems with large values near boundaries ]
   
   The problem was in the boundary condition routine, but also in advection. 
   In (2d) advection, remember to pass nx+1 staggered velocities to advect1d, and also remember to
   pass -2:nx+3 values of the scalar s (that is, interior + ghost points) to advect1d.
   

2. For pgm3, should we put the initial cone center to (0,0) or (0.0, -0.05)?
   Can [I] use the knowledge that the number of ghost points is 3 when I organize loops?
   
   I think this is consistent on the web page - finally!  Please use (0.0, -0.05).

3. Can I "hardcode" the number of ghost zones at the recommended value of 3?
   
   Yes, this is fine.

   ___________________________________ older questions/answers _________________________________________
   
4. I am getting a floating point exception ...
   
   Presumably this could be a divide-by-zero error but it also could be that your solution is blowing up.
   I recommend the following:

   a) run your code, plotting the output fields every time step.  Check the last few plots you were able to get, for any sign of the solution
       blowing up (this is more likely to occur near domain boundaries than in the interior)
   b) recompile your code, using compile options recommended on page 2 of the computing handout from the first class,
        under "to debug your program", and run again to see if any problems show up.
   c) print statements - or use of a debugger - are also useful.  I print out the min/max values in my code on every time step -
       you may find this helpful as well.

5. [Fortran] The time shown on the "sfc" plots becomes asterisks ***** on/after T=100 seconds
   
   If so, alter the format statement for writing the time, e.g.
   
         write(tlabel,15) time
   15  format('TIME =',f7.4)
   
   Change the format from f7.4 to say f9.4, to handle the larger time values. 
   I have made this change in my 502/Pgm3 copy of the sfc.f90 routine.
   
   
6. Should the Takacs scheme also have a symmetric solution (about X=0)?  Mine is symmetric for other schemes
   
   Yes it should.
   
   
7. [beware of checks for 'solution blowing up']
   
   A holdover from program 1 is checks for "unreasonable" values in your solution.  You may be checking to see if smin is below -5, or smax greater than 15.  These checks can be commented out or deleted here (plotting smin and smax is not required.  I still plot it in my code).  The solution - particularly when using Lax-Wendroff - can easily exceed these earlier thresholds.  If you want to check for a failing solution, be less restrictive, e.g. "if (smin.lt.-25.0 .or.  smax.gt.35.0) then..."
   
   
8. [more on codes failing with segmentation faults or behaving different from one run to another]
   
   Another way in which problems are occurring in Fortran is when updating your 2d array with results from the 1-d array returned by the 1-d advection routine.  Your 1-d advection routine only updates the 1d array from 1 to nx, e.g. "s1d_new(1:nx)".  If you attempt to copy boundary points - which were not updated by the 1d routine - from this 1d array, this could lead to errors or unpredictable behavior.  That is, "s1d_new(nx+1)" was not set by your 1d routine, and so should not be used.
   
   Explanation: this problem is most severe when the 1-d array in question is defined internally in a subroutine.  For example, say you pass the dimension "nx" to routine "advection".  Inside "advection" you dimension/allocate an array, "real, dimension (-1:nx+2) :: s1d_new", and then proceed to use s1d_new in calls to your 1d advection routine.  All well and good, but when you call advection, "s1d_new" is essentially allocated at that time.  Its contents, unless you explicitly set them, are undefined.  So when s1d_new is returned from advect1d, only the contents from 1 to nx have been set.  The other values could contain anything, and trying to access values outside [1,nx] may produce an error - or may not, depending on what was in memory when the program was running.  Some computers/compilers will set things to zero for you.  You shouldn't count on this.  Just because an array has been allocated does not mean its values are defined - unless you define them.
   
   
9. My code is blowing up unexpectedly and/or behaving differently when I run it different times
   
   Compilers can be unforgiving if you move or manipulate data you have not set previously.  For example, if you only set values for s1d(1:nx), and try to access s1d(nx+1) -- even if it is dimensioned OK - you may get errors or inconsistent behavior from one run to the next.  Here are suggestions that have helped some:
   
         a) compile with "-zero" (this refers to Fortran only).
             This initializes all values to zero (integers, real, everything) at the start of the program.
   
         b) set the entire contents of the s1() and s2() arrays to zero before calling your ic() routine. 
             In Fortran, you can just say "s1=0.0" and "s2=0.0".
   
   
10. I'm having (code) problems...
    
    * Fortran programmers:  Do take advantage of Fortran's additional debugging options. 
       Compile your code with "-check  all   -traceback"; it can find common problems with array dimensions or loop bounds. 
    * In C, use "-g  -debug  extended  -traceback".
    * For either: should you not get useful information (particularly, in what routine the program is having e.g. "segmentation fault"),
       try running your code within a debugger.  One way: use

          gdb  name-of-executable-program

       to run it; then type run, then proceed as usual. Type quit to exit the debugger.

11. Is the formula for the initial scalar cone "s" the same as in program 2?
    
    Yes, use the same formula.  Note the location of the center, and the cone radius, differ from program 2.
    
    
12. Won't the physical location of U and V now extend beyond the [-0.5,+0.5] range of the domain we are using?
    
    Yes.  The range of X and Y between [-0.5,+0.5] applies only to the scalar field S.  Since U is staggered 1/2 grid length to the left of S, and since there are nx+1 values of U vs. nx values of S, U(1) will be dx/2 to the left of S(1) and U(nx+1) will be dx/2 to the right of S(nx).  Similarly, V (from j=1 to ny+1) will also extend outside the [-0.5,+0.5] physical range.
    
    
13. Say more about the x, y location differences for s, u, and v ...
    
    In Program 2, we could get away with using the same x,y locations for the scalar and flow field because U was a function of only Y, and V was only a function of X.  This is no longer true.  In fact, you should now be using three separate loops for your initial conditions for S, U, and V, because the array sizes and thus do- or for-loop ranges differ between S, U, and V.  Relative to the X, Y location of a S point, the location of a U point is (X-dx/2, same Y); the location of a V point is (same X; Y-dy/2).
    
    
14. The local Courant number calculations ...
    
    Since U and V are now functions of X and Y, we need to take into account the grid staggering. When computing the local Courant number for a S point for Lax-Wendroff, 6^th-order Crowley, or Takacs, average U in the x-direction, and V in the y-direction
      [remember, U(i,j) is dx/2 to the left of S(i,j), so the average U at a S point location is 0.5*(u(i,j)+u(i+1,j)) ]
    
    
15. I noticed that you have plots of wind vectors in your program 3 example figures. Would I be able to get access to this routine?
    
    Yes - this is via a NCAR Graphics routine, the simplest of which I am using here (it is not required).  I called it as:
    
    
    call ezvec(uplot,vplot,nx,ny)
    
    
    where uplot, vplot are my U and V wind components that have been averaged to the S locations.
    Documentation on ezvec is available here:  http://linux.die.net/man/3/ncl_ezvec
    Note: this routine is obsolete - has been replaced by others, documented here:  http://ngwww.ucar.edu/supplements/vectors/
    
    
    
16. My plots aren't looking like the test results..
    
    First compare the min/max stats, and note I used a contour interval of 2.0 here, instead of 0.5 from program #2. 
    I have added a note above to clarify this.