#
# scheme:
#
# A scheme is an algorithm for advancing the field in time.
# It consists of a stencil operator, left and right boundary conditions,
# and a procedure that returns the maximum time step allowed by stability.
#

printf("scheme procedures:  \
CIR_stencil, \
CIR_l_bc, \
CIR_r_bc, \
CIR_max_time_step, \n\
LF_stencil, \
LF_l_bc, \
LF_r_bc, \
LF_max_time_step, \n\
LW_stencil, \
LW_l_bc, \
LW_r_bc, \
LW_max_time_step, \n\
Godunov_stencil, \
Godunov_l_bc, \
Godunov_r_bc, \
Godunov_max_time_step, \n\
Osher_stencil, \
Osher_l_bc, \
Osher_r_bc, \
Osher_max_time_step\n");

printf("scheme global variables:  \
CFL_factor, \
artificial_viscosity, \
Riemann_solver, \
Osher_flux, \
CIR_scheme, \
LF_scheme, \
LW_scheme, \n\
Godunov_scheme, \
Osher_scheme\n");

#
# global variables
#

CFL_factor := .9:
artificial_viscosity := .1:
Riemann_solver := Burgers_Riemann_solver:
Osher_flux := Burgers_Osher_flux:
CIR_scheme := [CIR_stencil, CIR_l_bc, CIR_r_bc, CIR_max_time_step]:
LF_scheme := [LF_stencil, LF_l_bc, LF_r_bc, LF_max_time_step]:
LW_scheme := [LW_stencil, LW_l_bc, LW_r_bc, LW_max_time_step]:
Godunov_scheme := [Godunov_stencil, Godunov_l_bc, Godunov_r_bc, Godunov_max_time_step]:
Osher_scheme := [Osher_stencil, Osher_l_bc, Osher_r_bc, Osher_max_time_step]:


CIR_stencil := proc(u_l, u_m, u_r, dx, dt)

#
# Courant-Isaacson-Rees stencil operator
#

  local lambda;

  lambda := (dt/dx)*velocity(u_m);
  if (lambda > 0) then
    lambda*u_l + (1.-lambda)*u_m
  else
    (1.+lambda)*u_m - lambda*u_r
  fi
end:


CIR_l_bc := proc(u_m, u_r, dx, dt)

#
# Left boundary condition scheme for the Courant-Isaacson-Rees scheme.
#

  CIR_stencil(u_m, u_m, u_r, dx, dt)
end:


CIR_r_bc := proc(u_l, u_m, dx, dt)

#
# Right boundary condition scheme for the Courant-Isaacson-Rees scheme.
#

  CIR_stencil(u_l, u_m, u_m, dx, dt)
end:


CIR_max_time_step := proc(field, grid)

#
# Maximum time step for the Courant-Isaacson-Rees scheme.
#

  local max_c;
  global CFL_factor;

  max_c := max(op(convert(map(abs,map(velocity,field)),list)));
  CFL_factor*get_spacing_of_grid(grid)/max_c
end:


LF_stencil := proc(u_l, u_m, u_r, dx, dt)

#
# Lax-Friedrichs stencil operator
#

  .5*((u_l + u_r) - (dt/dx)*(flux(u_r) - flux(u_l)))
end:


LF_l_bc := proc(u_m, u_r, dx, dt)

#
# Left boundary condition scheme for the Lax-Friedrichs scheme.
#

  LF_stencil(u_m, u_m, u_r, dx, dt)
end:


LF_r_bc := proc(u_l, u_m, dx, dt)

#
# Right boundary condition scheme for the Lax-Friedrichs scheme.
#

  LF_stencil(u_l, u_m, u_m, dx, dt)
end:


LF_max_time_step := proc(field, grid)

#
# Maximum time step for the Lax-Friedrichs scheme.
#

  local max_c;
  global CFL_factor;

  max_c := max(op(convert(map(abs,map(velocity,field)),list)));
  CFL_factor*get_spacing_of_grid(grid)/max_c
end:


LW_stencil := proc(u_l, u_m, u_r, dx, dt)

#
# Lax-Wendroff stencil operator
#

  local dtdx, f_m, u_ml, u_mr, av_term;
  global artificial_viscosity;

  dtdx := dt/dx;
  f_m := flux(u_m);
  u_ml := .5*((u_l + u_m) - dtdx*(f_m - flux(u_l)));
  u_mr := .5*((u_m + u_r) - dtdx*(flux(u_r) - f_m));
  av_term := artificial_viscosity*(u_r - 2.*u_m + u_l);
  u_m - dtdx*(flux(u_mr) - flux(u_ml) - av_term)
end:


LW_l_bc := proc(u_m, u_r, dx, dt)

#
# Left boundary condition scheme for the Lax-Wendroff scheme.
#

  LW_stencil(u_m, u_m, u_r, dx, dt)
end:


LW_r_bc := proc(u_l, u_m, dx, dt)

#
# Right boundary condition scheme for the Lax-Wendroff scheme.
#

  LW_stencil(u_l, u_m, u_m, dx, dt)
end:


LW_max_time_step := proc(field, grid)

#
# Maximum time step for the Lax-Wendroff scheme.
#

  local max_c;
  global CFL_factor;

  max_c := max(op(convert(map(abs,map(velocity,field)),list)));
  CFL_factor*get_spacing_of_grid(grid)/max_c
end:


Burgers_Riemann_solver := proc(u_l, u_r)

#
# zero-speed Riemann solver for Burgers' equation
#

  local s;

  if (u_l > u_r) then
    # shock case
    s := .5*(u_l + u_r);
    if (s > 0.) then
      u_l
    else
      u_r
    fi
  else
    # rarefaction case
    if (u_l > 0.) then
      u_l
    elif (u_r < 0.) then
      u_r
    else
      0.
    fi
  fi
end:


traffic_Riemann_solver := proc(u_l, u_r)

#
# zero-speed Riemann solver for the traffic dynamics
# equation u_t + (u(1-u))_x = 0
#

  local s;

  if (u_l < u_r) then
    # shock case
    s := 1. - (u_l + u_r);
    if (s > 0.) then
      u_l
    else
      u_r
    fi
  else
    # rarefaction case
    if (u_l < .5) then
      u_l
    elif (u_r > .5) then
      u_r
    else
      .5
    fi
  fi
end:


Godunov_stencil := proc(u_l, u_m, u_r, dx, dt)

#
# Godunov stencil operator
#

  local f_l, f_r;
  global Riemann_solver;

  f_l := flux(Riemann_solver(u_l, u_m));
  f_r := flux(Riemann_solver(u_m, u_r));
  u_m - (dt/dx)*(f_r - f_l)
end:


Godunov_l_bc := proc(u_m, u_r, dx, dt)

#
# Left boundary condition scheme for the Godunov scheme.
#

  Godunov_stencil(u_m, u_m, u_r, dx, dt)
end:


Godunov_r_bc := proc(u_l, u_m, dx, dt)

#
# Right boundary condition scheme for the Godunov scheme.
#

  Godunov_stencil(u_l, u_m, u_m, dx, dt)
end:


Godunov_max_time_step := proc(field, grid)

#
# Maximum time step for the Godunov scheme.
#

  local max_c;
  global CFL_factor;

  max_c := max(op(convert(map(abs,map(velocity,field)),list)));
  CFL_factor*get_spacing_of_grid(grid)/max_c
end:


Burgers_Osher_flux := proc(u_l, u_r)

#
# Osher's flux for Burgers' equation
#

  if (u_l > 0. and u_r > 0.) then
    .5*u_l^2
  elif (u_l < 0. and u_r < 0.) then
    .5*u_r^2
  elif (u_l >= 0. and 0. >= u_r) then
    # Osher's entropy fix for a transonic shock
    .5*u_l^2 + .5*u_r^2
  else
    0
  fi
end:


traffic_Osher_flux := proc(u_l, u_r)

#
# Osher's flux for the traffic dynamics equation u_t + (u(1-u))_x = 0
#

  local s;

  if (u_l < u_r) then
    if (u_l < .5 and .5 < u_r) then
      # Osher's entropy fix for a transonic shock
      RETURN(u_l*(1.-u_l) + u_r*(1.-u_r) - .25)
    fi;
    # shock case
    s := 1. - (u_l + u_r);
    if (s > 0.) then
      u_l*(1.-u_l)
    else
      u_r*(1.-u_r)
    fi
  else
    # rarefaction case
    if (u_l < .5) then
      u_l*(1.-u_l)
    elif (u_r > .5) then
      u_r*(1.-u_r)
    else
      .25
    fi
  fi
end:


Osher_stencil := proc(u_l, u_m, u_r, dx, dt)

#
# Osher stencil operator
#

  local f_l, f_r;
  global Osher_flux;

  f_l := Osher_flux(u_l, u_m);
  f_r := Osher_flux(u_m, u_r);
  u_m - (dt/dx)*(f_r - f_l)
end:


Osher_l_bc := proc(u_m, u_r, dx, dt)

#
# Left boundary condition scheme for the Osher scheme.
#

  Osher_stencil(u_m, u_m, u_r, dx, dt)
end:


Osher_r_bc := proc(u_l, u_m, dx, dt)

#
# Right boundary condition scheme for the Osher scheme.
#

  Osher_stencil(u_l, u_m, u_m, dx, dt)
end:


Osher_max_time_step := proc(field, grid)

#
# Maximum time step for the Osher scheme.
#

  local max_c;
  global CFL_factor;

  max_c := max(op(convert(map(abs,map(velocity,field)),list)));
  CFL_factor*get_spacing_of_grid(grid)/max_c
end: