###########################################################
# Wed Jun 25 18:39:08 PDT 2003
#
# This code solves the 1D Burger's equation
# using a Roe solver. See
#
#  http://laplace.physics.ubc.ca/pi-nr/lec/week2/Lec4.pdf
#  
# starting at page 10 for documentation.  
###########################################################

# Set the memory size (only needed for Fortran)
system parameter int memsiz := 12000000

# Domain extrema
parameter float xmin := 0.0
parameter float xmax := 5.0

# Number of Ghost cells used. 
parameter  int   Ng       := 2

# Parameters to define the initial conditions.

parameter int   ini_type  := 1
parameter float q_xc      := 2.5e0
parameter float q_R       := 0.5e0
parameter float q_L       := 0.5e0
parameter float q_0       := 2.5e0
parameter float q_amp     := 1.0e0
parameter float q_delta   := 0.5e0

# Define coordinate system:

cartesian coordinates t,x

uniform cartesian grid g1 [1:Nx] {xmin:xmax}

#Declare grid functions:

float q            on g1 at 0,1 { out_gf := 1 }
float q_nc         on g1 at 0,1 { out_gf := 1 }
float q_c          on g1 at 0,1 { out_gf := 1 }
float q_2c         on g1 at 0,1 { out_gf := 1 }
float FF           on g1 at 0   { out_gf := 0 }
float res          on g1 at 0   { out_gf := 1 }

############################################################
# Operator
############################################################
operator D_CN(f,t)     := (<1>f[0] - <0>f[0])  / dt
operator D_UPWIND(f,x) := (<0>f[0] - <0>f[-1]) / dx
operator MU(f,t)       := (<1>f[0] + <0>f[0] ) / 2.0
operator D_BW(f,x)     := ( 3 * <0>f[0] - 4 * <0>f[-1] + <0>f[-2])/(2 * dx)
############################################################
# Residuals
############################################################

evaluate residual q_nc {
   [1:1]   :=  D_CN(q_nc,t) = 0;
  [2:Nx]   :=  D_CN(q_nc,t) + MU(q_nc*D_UPWIND(q_nc,x),t);
}

evaluate residual q_c {
   [1:1]   :=  D_CN(q_c,t) = 0;
  [2:Nx]   :=  D_CN(q_c,t) + MU(D_UPWIND(0.5*q_c*q_c,x),t);
}
evaluate residual q_2c {
   [1:2]   :=  D_CN(q_2c,t) = 0;
  [3:Nx]   :=  D_CN(q_2c,t) + MU(D_BW(0.5*q_2c*q_2c,x),t);
}
#looper standard
looper iterative
auto update q_nc, q_c, q_2c

step.inc step updates q, res 
   header q, FF, lambda, x, Ng, res, t


init_q.inc init_q initializes q
   header q, q_xc, q_R, q_L, q_0, q_amp, q_delta, x, ini_type

initialize q_nc {
   [1:Nx]:= q
}  

initialize q_c {
   [1:Nx]:= q
}  

initialize q_2c {
   [1:Nx]:= q
}  

auto INITIALIZE q_nc, q_c, q_2c