Description Problem:  nonlinear quasi-compressible flow in 2D Due: 8 pm Wednesday April 17 -- entirely online on Moodle (nothing printed) Methods: Spatial advection:  Forward-time Lax-Wendroff and Piecewise-linear (scalar), box (motion variables) Diffusion: Leapfrog-time (momentum), time-lagged 2nd-order.  Pressure and buoyancy: Leapfrog-time, forward-backward. Domain:  2-D (vertical slab) Problem handout (for programs 5+6): PDF Hand-in settings for Program 5 Configuration: 451x51 domain (nx,nz); grid spacing 50 meters in x, z Run:  to 3500 time steps (700s physical time) using time step dt = 0.20 sec. Initial theta perturbations: (each are cold, -25 degree perturbation superimposed on 300K background theta value) Perturbation #1:  centered at x=25, z=1025 meters (centered on the left edge) Perturbation #2:  centered at x=12025, z=1025 meters  (not all the way to the right edge) Both perturbations have "radii" (denominators) of 2500 m (horizontal), 750 m (vertical) Advection:  Lax-Wendroff only, not piecewise-linear. Diffusion: coefficient for U and W = 20; coefficient = 45 for theta. Lax-Wendroff solution is noisier than piecewise-linear would have been...with PL we could have used less theta diffusion. Continuity equation: reduced sound speed Cs = 80 meters/sec Plotting: plot Theta' (the perturbation theta-thetabar) every 500 steps (every 100s physical time) with a contour interval of 1 degree C plot U, W and P only at 1500 and 3000 steps (300s and 600s physical time) with contour intervals:  2 (for U), 2 (W), 25 (P'). No hardcopy needed at all. Upload plots (as images e.g. GIFs or PNGs) and code (raw code files + makefile) on Moodle. If your code is blowing up, try higher diffusion to see if it will work then (and let me know what you used).  If need be, decrease dt but adjust how often you plot so your plots are near the physical times I have requested. If your code refuses to run all the way to 3500 steps, plot all fields at least every 500 steps using the original code settings above. Suggested order of development:  Computer Program #5 Start with a copy of program #2. Change the code to permit NX not equal to NZ. Change physical dimensions, so domain no longer [-0.5 : +0.5] Create initial theta field; plot (theta-<theta-bar>); verify looks OK e.g. compare to initial theta plot in Test A (PGF-only test), below. Dimension all arrays you will need (not including ghost points here) u, w, p: three time levels, 2-D (nx,nz; u has 1 extra point in X, and w has one extra in Z) theta: two time levels, 2-D (nx,nz) 1-D base-state density (and density @ w levels): 1-D (nz) Boundary condition routine Carry out pressure gradient force test (Test A) Carry out linear advection tests for theta (Test C, and test B if you want to check your boundaries) Carry out diffusion test for theta Full physics (Test D) If your Test D may have a U/W advection problem, try test E. Test case results for Program #5 Use these test cases to check your code at dx,dz = 100 meters and for advection, diffusion, and pressure/buoyancy processes. All test cases are ready via links below. Click on the small image to jump to a test page. Test A PGF / buoyancy test Test B T advection: splitU,W Test C T advection: U,W const. Test D Full physics, sinking thermal Test E [optional] Nonlinear advection test Update step: what needs to be done Your update step - after your integration (advection, diffusion, pressure terms) are done - will look like - If (this is the end of the very first time step) Then     Set tstep = 2*dt                               (it was previously = ?t before the time step loop started)     u2=u3, w2=w3, p2=p3, t1=t2       (copies n+1 time level over n; n-1 level unchanged) Else     u1=u2, w1=w2, p1=p2                 (copies n time level over n-1 level)     u2=u3, w2=w3, p2=p3, t1=t2       (copies n+1 time level over n) Endif Reference information: Array dimensions, indices, and boundary conditions  Details follow.  You really need only refer to this reference as questions arise. Overview: I dimension u, w, theta and p as (0:nx+1,0:nz+1), but have one additional point for each direction in theta in my 1-d passes directly.  Theta needs (-1:nx+2, -1:nz+2) if you dimension your theta array to include ghost points for the Piecewise Linear scheme.  P doesn't really need any ghost points (but I dimension it with one ghost point anyway). What array dimensions do you need?  Looking first at variables on the staggered grid... U (horizontal flow, meters/sec) is integrated from 2:nx, 1:nz; values at i=1 and i=nx+1 are set from symmetry BCs W (vertical flow, meters/sec) is integrated from 1:nx, 2:nz; values at k=1 and k=nz+1 are fixed at 0 for all time Theta (potential temperature, degrees Kelvin) is integrated from 1:nx, 1:nz Pressure (perturbation, p', in Pa) is integrated from 1:nx, 1:nz Physical processes Advection: U advection requires U values at i-1, i, i+1, and k-1, k, k+1 W advection also requires W values at i-1, i, i+1, and k-1, k, k+1 Remember U and W advection is done directly - no directional splitting is used. Theta advection for Piecewise Linear requires values (in your 1-D passes) from -1:nx+2 (Fortran notation).  Why? Your solution for i=1,...,nx requires fluxes at (i-1/2) and (i+1/2), which we reference as F(i) and F(i+1) Fluxes at i-1/2 (i.e. F(i)) require the scalar (theta) at i and i-1, and slopes at i and i-1. Thus you would need theta from i=0 (i-1 point used by F(i)) to i=nx+1 (i point used by F(i+1)) But slopes at i and i-1 require scalars one point on either side, so - you need theta from i-2 (for F(i), slope(i-1) left point) to i+2 (for F(i+1), slope(i) right point) This last criteria means you need 2 ghost points for theta, from -1:nx+2 and -1:nz+2. I dimension my theta arrays 0:nx+1,0:nz+1 like the rest, but handle the one extra point (needed for the 1-d piecewise linear passes) in the 1-d advection routine directly.  It would be better to dimension the theta arrays as -1:nx+2 and -1:nz+2, but my 'shortcut' is acceptable in this case. Diffusion: For U, W, and theta: simple 2nd-order derivative requires i-1,i,i+1 and k-1,k,k+1 I did diffusion everywhere by which I mean my code contains: mixing on U from k=1 to k=nz (but still from i=2,nx: you set i=1 and i=nx+1 from BCs) mixing on Theta from k=1 to k=nz mixing on W from i=1 to nx but only k=2 to k=nz, since W(k=1) = W(k=nz+1) = 0 always. Pressure terms, and buoyancy: For pressure gradient term for U, you need dp/dx.  So u(i,k) requires p(i,k) and p(i-1,k). For pressure gradient term for W, you need dp/dz.  So w(i,k) requires p(i,k) and p(i,k-1) ... as well as density averaged to the same W height as W itself. The buoyancy term for W requires an average of theta' (theta-thetabar) at levels above and below the W height; for w(i,k), you therefore need theta at k and k-1. For the pressure equation itself, you need: density(k), and du/dx:  du/dx for a p location requires u at i+1 and i. density at W levels, and dw/dz - denoting density at W levels as density_W: d/dz( density*W ) centered on a p location requires density_W and W at k+1 and at k. Density terms: You only need 1-D arrays (dimensioned with height) for density. It is not evolved with time.  "rho-overbar" is a function of height only. It is simplest to compute density averaged to W levels ahead of time, and carry around two 1-D arrays, density(nz) and densityW(nz+1).  Note that you can't average density to W height k=1 and k=nz+1, since there is no density below k=1 or above k=nz.  But these values (density_W at k=1 and at k=nz+1) are always multiplied by W at that same level, and W is zero at the very top and bottom.  So -- I set my averaged density_W at k=1 and at k=nz+1 to 9999 ... it shouldn't matter; it will be multiplied by zero. Boundary conditions This is covered in the handout in detail.  But as a summary: Theta:  you need theta at i=0 and at i=nx+1 and at k=0 and k=nz+1 for diffusion. You need an extra theta ghost point (as noted above) for 1-d advection passes with the PL routine. U: in the horizontal you only evolve 2:nx, and set i=1 and i=nx+1 from BCs. But in the vertical you'll need a ghost point for diffusion and for the advection of U. W: you evolve 2:nz, and set k=1 and k=nz+1 to zero for all time.  You'll need ghost points in the horizontal for diffusion and for advection of W. P': no ghost points really needed, given the way we handle U and W. X symmetry boundaries: u changes sign; theta, p and w keep same sign across the boundaries. Z zero-gradient boundaries: all variables keep same values above/below the physical domain as at their boundaries

Frequently asked Questions & Answers, and Hints

1. Did you generate the test cases using Fortran?  I am using C, the array index starts from 0.  My plots look exactly the same
   except my x and z coordinates are one [larger-?] than yours.  For example, my theta plot at initial condition says the min is at (100, 20) and the max is at (0, 0) but your graph says the min is at (101, 21) and the max is at (1, 1).  These differences would probably make sense if you generate[d] them using Fortran.
   
   All test cases were from Fortran - typically this means you would subtract one from each index to match C results.  But the coordinate depends on the number of ghost points you are using, and the fact that Fortran arrays start at 1 while C starts at 0.   You might check primarily to see that the theta temperature perturbation is centered in your domain, based on the indices you are using.
   
   
2. What is the difference between PL monotonic and non-monotonic?  As far as the description on program3 goes, it says what makes
   the function non-monotonic is deltaq.  But what do we need to make it monotonic?
   
   The only difference is the monotonic filter we discussed in class - where slopes (the scalar centered derivative dq/dx) are set to zero
   if there is a minimum or maximum in the zone.  But note you are asked to use the non-monotonic version for programs 5 and 6.
   
   
3. On page 7 of the handout, for the BC routine, it implies that we need 0-gradient boundary conditions for U for the top and bottom boundaries.
   I thought the size of U is (1:nx+1, 1:nz).  I may need to increase the size to include ghost points for U?
   
   Yes, U, W and P all need ghost points.  So include at least 1 ghost point above and below the physical points.
   
   
4. [Also] page 7 of the handout implies we need to make W same on either side of the symmetry boundary.  I didn't get this because W lies along the symmetry boundary (dashed yellow lines).  I thought the size of W is (1:nx, 1:nz+1).  I need to add ghost points for W?
   
   Yes, W needs ghost points.  If you are evolving W on the left edge (i=1 in Fortran), and you need W at i-1, you will use symmetry boundary
   conditions set at the ghost point for (again, Fortran) i=0.  Our symmetry BCs for W say that W at i=0 is equal to W at i=2.
   
   
5. The perturbation pressure equation should be used inside the PGF function, although it is not highlighted in slide 10 of class #25, right?
   
   Right!  My mistake.  For the PGF routine (really pressure gradient+buoyancy routine):
   
   ? The PGF routine should have passed to it u3, w3, p1, p3, and t1 [some may call this q1]. 
     That is, we pass the "latest" arrays (the n+1 time level arrays) for u, w, and p, along with the n time level of potential temperature
     (t1 or q1, whatever you call it), and the prior (n-1) time level of p. 
   
   ? The order of tasks in the PGF routine is: update u; update w, set BCs for the updated u and w, and then update p.
   
   ? For the flow field, we set u3 = u3 + (pressure gradient terms), and w3 = w3 + (pressure gradient terms) + (buoyancy terms).
      Remember to use the (n-1) time level of pressure "p" when computing the pressure gradient terms for u, w.
   
   ? For the pressure field, you compute du/dx and d(density*w)/dz using your just-updated (n+1)-time level arrays for u and w.
   
   
6. Why do we need three "p" arrays?  In PGF we solve p3=p1+.... can't we just use one (or at most two) arrays for pressure?
   
   We really do need three pressure arrays.  The centered time differencing has us jumping from (n-1) to (n+1), but the update step after our pressure gradient routine requires copying (n) to (n-1), and then (n+1) to (n).  We don't use the n-time-level pressure array anywhere, but we do shuffle values among those three time levels.
   
   
7. Likewise, why do we now need two potential temperature arrays?  I've used just one in programs 2-4.
   
   We haven't really needed more than one scalar ("q" in programs 2-4; I called it "t" in program 5, you choose) array until now.  The difference now is that after we finish our advection routine, we still need the old (time level "n") potential temperature for use in the diffusion routine (for diffusion of theta) and in the buoyancy term for "w" in the PGF routine.  So advection can't simply overwrite one potential temperature array.  So as you code things, keep the n-time-level (current) array for theta unchanged, for use in diffusion and PGF, and use your (n+1)-time-level array for theta separate as you update it with the advection and diffusion contributions.
   
   
8. My (U,W, and/or P') fields are remaining at zero instead of evolving.
   
   The first place to look is probably your update step.  Is it possible you are copying an array of zeroes over your other array time levels?  If not that, note that if you have a statement like:
   
       u3 = u1 - ....
   
   in your advect routine, for example, and you comment out the call to advect to do a PGF test, you will need to add the statement "u3=u1" before calling PGF.  Otherwise, arrays like U and W may not be correctly taking the leapfrog step forward.  As a matter of habit, I put u3=u1  and  w3=w1   in my main program before calling advect, diffusion, and pgf, and those subroutines then just add to the new (n+1) values; they have statements like u3 = u3 + ...  and  w3 = w3 + ... in the advect, diffusion, and PGF routines.
   
   
9. Do we need to take 2 theta time steps for every leapfrog U/W/P' time step?
   
   No.  We can 'mix' the forward-in-time theta integration and leapfrog-in-time U/W/P integration without taking double the number of theta steps.  Remember the leapfrog time is jumping 2*dt, but is doing so from (n-1) to (n+1), whereas theta is going from (n) to (n+1).
   
   
10. On page 1, for the discrete form of equation for 'theta', there are two versions given. Why the first one is given in matrix form, does the upper entry in the matrix (big bracket) means in X-direction,  lower entry means Z-direction?
    
    Yes, exactly.  Really, other than boundary conditions, you don't change your scalar (theta) advection at all for this program.  The diffusion is off in a separate subroutine (with leapfrog rather than forward time differencing). 
    
    
11. Initial state: it says 'set rho_w level @ k=1 to k=9999 as a check' on handout p. 3, what do you mean by that? then set @ what value?
    
    What I meant is:  we create the density array I call rho at the U,Theta,P levels from the equations in the handout.  Then we average this density to the intermediate vertical velocity levels - I call this averaged density array "rho_w".  The very first rho_w (at k=1) is used in conjunction with vertical velocity, which is always zero at k=1.  So it makes sense to just set rho_w(k=1) = 9999 or some other large value, because it should always be multiplied by zero.  Otherwise it really isn't used, since W at k=1 is not evolved in time.
    
    
12. I have compared my results with yours and for time step 100 my u,w and p have similar results to yours. But, for time step 200 the u and w are a little bit different from yours but the pressure has values of 300-400 and your pressure field is about 600.  I dont know if this is a problem.
    
    I worry about pressure fields if the magnitude is increasing dramatically and consistently.  If your field is different from mine but the pattern is similar (after all, only the gradients matter), I wouldn't worry about it.  Only U, W, theta really matter - and the pattern of P should be similar (otherwise the other fields wouldn't be).
    
    So, no, I won't be grading on absolute values of P.  It would be incorrect as we have no real constraint on P', and I'm sure coding and compiler differences could lead to a drift.  Perhaps if everyone removed the mean P' on every time step, it might make sense to compare the values instead of just the pattern/gradients, but as it is, it doesn't make sense to do so.
    
    
13. [Numerous questions about early vertical velocity field magnitudes in the PGF test, test A]
    
    I believe my results are correct.  If you find your vertical velocity amplitudes are appreciably larger (say, the W field minimum at step #2) than mine, you might check to make sure your (1) array updates are being done correctly, (2) that you set u3=u1 and w3=w1 at the start of each time step (that way you can add in contributions from various processes with w3=w3+...  and  u3=u3+...    ; and (3) that your "tstep" (my time increment for u, w, and p') doubles from dt to 2*dt after the first time step, to transition your code from forward in time to leapfrog time for u, w, and p'.
    
    
14. [Boundary problems in scalar field in the advection tests]
    
    The problem turned out to be with boundary conditions in the scalar advection code; theta BCs were not being set between the directionally split advection steps (between X and Z, or Z and X).  In their case the problem did not occur if the mean potential temperature (theta-bar) was 0, but was apparent if it was = 300K.
    
    
15. [Looking at text output from test case B]: "you set dt to 2dt after initial; but each timestep represents dt.  Why are we increasing dt to 2dt again?"
    
    That file listing includes the following:
    
    

      step    time     umin      umax      wmin       wmax     tmin     tmax       pmin       pmax
    ------------------------------------------------------------------------------------------------
        0    0.000 -20.00000  20.00000   0.00000   0.00000   0.00000  25.00000    0.0000   0.0000
     Changing to 2dt =  0.5000000    
        1    0.250 -20.00000  20.00000   0.00000   0.00000  -0.00092  24.99994    0.0000   0.0000
        2    0.500 -20.00000  20.00000   0.00000   0.00000  -0.00174  24.99966    0.0000   0.0000

    The change shown was not dt, but tstep.  Tstep was changed from dt, for the first step, to 2*dt, which is used
    for all other time steps in the u,w,p expressions which look like:
    
             u3 = u3 + tstep*(terms)
    
    
16. [Questions about staggering and averaging, in the pressure gradient terms; do we need to average p' since not at u,w locations?]
    
    True, the p locations are different from u and w, but note the terms are really:
    
    u:  dp'/dx
    w:  dp'/dz
    
    both of which, due to the staggering, are actually valid at u and w locations; dp/dx is centered on U, and dp/dz is centered on W.
    No additional averaging is necessary.  For the U and W equations, the pressure terms we use are derivatives valid at the U and W locations.
    
    
17. [Questions about page 5, PGF subroutine bullets a and b]
    
    There is a correction: both should have "tstep*(pgf terms)" in those expressions.
    
    
18. [Questions about theta diffusion]
    
    The governing equations at the top of page 1 do not (but should) list diffusion for theta, though they are in the discretized equations.
    The discretized equations are correct.
    NOTE there is theta diffusion; your diffusion will not vary depending on the advection form you use (the handout erroneously said "only if Lax-W").
    
    
19. Why are the units for the density g kg^-1? Why isn't it kg m^-3?
    
    You are right and the handout is incorrect.  It is kg/m**3.
    
    
20. What value of Cs (sound speed) we should use?
    
    I will specify it along with the hand-in problem settings.  For now, check against the test case results posted online (where it says what sound speed to use)
    
    
21. What is "rho_{theta level}" ??
    
    This is density (rho-bar) valid at the heights (above the lowest W level, where z=0) of the u, p and theta variables. 
    The height of W levels is dz/2 from those of u,p,theta.
    
    
22. I don't understand the top paragraph of page 3, where it says "Check the initial state with this data for dx=dz=100m, at k=11,..."?
    ... is this the initial condition we should use, with perturbations introduced only at level k=11? (plus the temperature perturbation
    introduced later on) ?
    
    No, that information is just for you to check your density calculations.  The temp perturbations are independent of this.
    
    
23. What is p' ? Is it P = overbar{p} + p' where overbar{p} is calculated using the eqn on page 2? and total P using the value from the first line
    on top of page 3? then we get the initial p'?
    
    We only use overbar P to calculate density.  Otherwise it is not needed.
    We set the initial p' field to zero, and let it evolve - we do not calculate or need the total pressure at all.
    
    
24. What are xradius and zradius?
    
    xradius and zradius are the terms in the denominator of the radius calculation expression (x and z terms).
    Poor choices of wording, probably, since radius implies a circular rather than elliptical perturbation for our initial scalar theta field.
    
    In programs 2/3, the initial condition expression used tests of d vs. r; here, the ratio of distance and radius is built into the r(m) expression, such that for r>1, there is no perturbation added to the initial theta (initial value = theta-bar = 300K) field, and for r <= 1, we use the expression in the handout.
    
    
25. What is directionally split, besides theta advection?
    
    
    Nothing.  Your perturbation temperature (theta) advection is the same as before -- directionally split.  Everything else -- theta diffusion; u, w diffusion; u, w advection; pressure gradient and buoyancy terms - are unsplit, evaluated directly in 2D.  So, e.g., the u advection looks like:

         u3 = u3 + tstep*[(x advection of u)+(z advection of u)]
    

26. What do you mean "free slip" boundary conditions for u ?
    
    
    Here it means you don't have to do anything.  The horizontal flow u "slips" freely above the ground surface located where w(k=1)=0.0 ... had we used "no slip" conditions for u, you might, for example, set u to zero at k=1 for all time.  The only relevant boundary conditions on u are how we set the ghost zones: 0-gradient at the top and bottom of the model domain.
    
    
    
27. [Not sure how to set up initial conditions for temperature with loop m=1,2 on handout page 3]
    
    
    For each location (i,k), you set t(i,k) to tbar, that is, to 300.0.  Then, include a loop from 1...2 that checks to see if your (i,k) grid point is within radius of 1.0; if so, add to your t(i,k) the appropriate contribution from temperature perturbation #1 or #2 using the cos(r) function on page 3 of the handout. 
    
    
    
28. [Hint: be very careful with array dimensions now that NX is not equal to NZ!]
    
    
    If you program in C, you may still have arrays dimensioned [NXDIM][NXDIM] or the like, and you may have for() loops ranging from i1...i2 for both "x" and "y" dimensions from your code in program #3.  So you will want to check carefully your array dimensions, and your do- (or for-) loop bounds in coding up program #5.
    
    
    
29. There are places in the handout, page 1, that show diffusion of theta, or of the perturbation theta' ... does it matter?  Is this an error?
    
    
    The handout is written in the more general case where the theta-bar value is a function of height: tbar(z), in my code.  In such a case, doing horizontal diffusion can be on theta or theta', since the mean will drop out (not affect) the calculation.  When tbar is a function of height, the mean will not "drop out" and so the diffusion we use is on the deviation from the mean, T-Tbar.  For that reason, I used d**(theta')/dz**2 in the handout, but since Tbar is not a function of height - but, rather, a constant - it doesn't matter whether your x or z diffusion is on the total theta (i.e. t1 array) or on the deviation (t1-tbar).
    
    
    
30. [Problems with boundary conditions for theta, as applied to advection of temperature]
    
    
    You can always enforce your temperature boundary conditions in the routine that calls your 1-d advection method.  If you are performing x-advection of theta, you can set the ghost zone values in the 1-d array q before you pass it to the 1-d advection (PL) routine.  Note that we have changed x-boundary conditions for theta since program #3; you may still have old 0-gradient code in there.  For vertical advection of theta, you can similarly set your ghost zone values of your 1-d q array to be 0-gradient before calling your 1-d advection routine.
    
    If not the above, then call a routine to set your temperature boundary conditions prior to calling your advection routine for theta.
    
    
    
31. Next?