#! /usr/bin/env
#
def lesson03 ( ):

#*****************************************************************************80
#
## lesson03 presents the issues of convergence and the CFL condition.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    14 November 2015
#
#  Author:
#
#    Initial version by Lorena Barba
#    This version by John Burkardt
#
#-------------------------------------------------------------------------------
#  Did you experiment in Steps 1 and 2 using different parameter choices? 
#  If you did, you probably ran into some unexpected behavior.  Did your 
#  solution ever blow up?  (In my experience, CFD students love to make 
#  things blow up.)  
#
#  You are probably wondering why changing the discretization parameters 
#  affects your solution in such a drastic way.  This notebook complements 
#  our interactive CFD lessons by discussing the CFL condition.
#  You can learn more by watching Professor Barba's YouTube lectures.
#
#  Convergence and the CFL Condition
#
#  For the first few steps, we've been using the same general initial and 
#  boundary conditions.  With the parameters we initially suggested, the grid 
#  has 41 points and the timestep is 0.25 seconds.  Now, we're going to 
#  experiment with increasing the size of our grid. 
#
#  The code below is identical to the code we used in Step 1, but here it has 
#  been bundled up in a function so that we can easily examine what happens 
#  as we adjust just one variable: the grid size.
#-------------------------------------------------------------------------------
#
#  numpy is a library for array operations akin to MATLAB
#
import numpy        
#
#  matplotlib is 2D plotting library
#         
from matplotlib import pyplot    
%matplotlib inline

def linearconv ( nx ):
#
#  NX is the number of nodes to use over the interval [0.0, 2.0].
#
#  DX will be the spacing between the nodes.
#
  dx = 2.0 / ( nx - 1 )
#
#  NT is the number of timesteps.
#
  nt = 20
#
#  DT is the size of each timestep.
#
  dt = .025
#
#  V is the constant velocity.
#
  v = 1.0
#
#  Set the X values of our nodes.
#
  x = numpy.linspace ( 0, 2, nx )
#
#  Define a numpy array of length NX, with every value set to 1.
#
  u = numpy.ones ( nx )
#
#  Values of U corresponding to 0.5 <= X <= 1 should be set to 2.
#  Python ranges stop just BEFORE the last element in the list, so
#  we have to write our range as ".5/dx : 1/dx+1" in order to pick
#  up the last element at 1/dx.
#
  u [ .5/dx : 1/dx+1 ] = 2.0
#
#  Make a second array, UN.  This will save the current value of U, while
#  we compute the next one.  We can set this to zero initially.
#
  un = numpy.zeros ( nx )
#
#  To start out, U contains our estimated solution at the 0th time.
#  Take NT time steps, and update the value of U each time.
#
  for n in range ( nt ):
    un = u.copy()
    for i in range ( 1, nx ):
      u[i] = un[i] - v * dt / dx * ( un[i] - un[i-1] )

  pyplot.plot ( x, u );

  return
#-----------------------------------------------------------------------------80 
#  End of the linearconv() function!  
#-----------------------------------------------------------------------------80
#  Now we can compute our numerical solution to the convection problem.
#  We assume that we can improve our results by using a finer mesh.
#  Note that, in this procedure, we are keeping the time step size fixed.
#-------------------------------------------------------------------------------
linearconv ( 41 )
#-------------------------------------------------------------------------------
#  This is the same result as our Step 1 calculation, reproduced here for reference.
#-------------------------------------------------------------------------------
linearconv ( 61 )
#-------------------------------------------------------------------------------
#  Here, there is still numerical diffusion present, but it is less severe.
#-------------------------------------------------------------------------------
linearconv ( 71 )
#-------------------------------------------------------------------------------
#  Here the same pattern is present -- the wave is more square than in the previous runs.
#-------------------------------------------------------------------------------
linearconv ( 85 )
#-------------------------------------------------------------------------------
#  This doesn't look anything like our original hat function.  What happened?
#
#  To answer that question, we have to think a little bit about what we're 
#  actually implementing in code.  In each iteration of our time loop, we use 
#  the existing data about our wave to estimate the speed of the wave in the 
#  subsequent time step. Initially, the increase in the number of grid points 
#  returned more accurate answers. There was less numerical diffusion and the 
#  square wave looked much more like a square wave than it did in our first 
#  example.
#
#  Each iteration of our time loop covers a time-step of length dt, which 
#  we have been defining as 0.025
#
#  During this iteration, we evaluate the speed of the wave at each of the x 
#  points we've created. In the last plot, something has clearly gone wrong.
#
#  What has happened is that over the time period dt, the wave is 
#  travelling a distance which is greater than dx. 
#
#  That is, the grid spacing has gotten so small that, in a single time step, 
#  the wave should travel more than the length of one grid cell.  But our
#  numerical scheme only allows information to travel at most one grid cell
#  width per step.
#
#  The length dx of each grid box is related to the number of total points nx, 
#  so stability can be enforced if the dt step size is calculated with respect 
#  to the size of dx in such a way that we have
#
#    v dt dx = c <= cmax
#
#  where v is the speed of the wave; c is called the CFL number (named for
#  Courant, Friedrichs, and Lewy) and the value of cmax that will ensure 
#  stability depends on the discretization used.
#
#  Another way of looking at this formula is to assume that V is fixed,
#  that we have chosen DX, and now must choose a time step DT that keeps
#  our C value below the limit:
#
#    dt <= cmax / v / dx
#
#-------------------------------------------------------------------------------
#  In a new version of our code, we'll use the CFL number to calculate the 
#  appropriate time-step dt depending on the size of dx.
#-------------------------------------------------------------------------------
import numpy
from matplotlib import pyplot

def linearconv2 ( nx ):
  dx = 2/(nx-1)
  nt = 20
  v = 1.0
  sigma = 0.5
#
#  Notice that now DT will vary
  dt = sigma * dx

  u = numpy.ones(nx) 
  u[.5/dx : 1/dx+1]=2

  un = numpy.ones(nx)

  for n in range(nt):  #iterate through time
    un = u.copy() ##copy the existing values of u into un
    for i in range(1,nx):
      u[i] = un[i]-v*dt/dx*(un[i]-un[i-1])
        
  pyplot.plot(numpy.linspace(0,2,nx),u)

  return

#-------------------------------------------------------------------------------
#  End of linearconv2 ( ).
#-------------------------------------------------------------------------------
#  Now let's look at how the revised code handles 20 time steps at a
#  sequence of finer and finer meshes.
#-------------------------------------------------------------------------------
linearconv2(41)

linearconv2(61)

linearconv2(81)

linearconv2(101)

linearconv2(121)
#-------------------------------------------------------------------------------
#  Notice that as the number of points nx increases, the wave convects a 
#  shorter and shorter distance. The number of time iterations we have advanced 
#  the solution at is held constant at nt = 20, but depending on the value of 
#  nx and the corresponding values of dx and dt, a shorter time window is being 
#  examined overall.
#
#  In a more typical computation, we would have a target time we wanted to reach.
#  In that case, once we have computed the appropriate DT, we would need to 
#  compute the value of NT, the number of timesteps, necessary in order for us
#  to reach the time T.
#
#-------------------------------------------------------------------------------
#  LEARN MORE:
#
#  It's possible to do rigorous analysis of the stability of numerical schemes.
#
#  Watch Prof. Barba's presentation of this topic in Video Lecture 9 on You Tube.
#-------------------------------------------------------------------------------