############################################################
# Solves Carlos Palenzuela et al's model for "stiff"
# wave/advection equation
#
#    w_t = v_x                                (1)
#
#    v_t = w_x + (-v + c w) / epsilon         (2)
# 
# on 0 <= x <= 1; t >= 0 with boundary conditions 
#
#    v(t,0) = v(t,1) = 0
# 
# and initial conditions 
#
#    v(0,x) = amp exp(-(x - x0)^2/delta^2)
#
#    w(0,x) = (v + epsilon exp(-(x - 0.5)^2/0.01) / c
#
# so that for epsilon -> 0, initial conditions are consistent 
# with pure left mover, while for epsilon bounded away from 
# 0, solution has both left and right movers
#
#  Difference scheme 
# 
#    o Crank-Nicholson with standard second-order spatial 
#      discretization, except at boundaries where O(h)
#      forward/backwards differences are used for v_x (2)
############################################################

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

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

# Gaussian initial data parameters 
parameter float amp        := 1.0
parameter float x0         := 0.5
parameter float delta      := 0.1

parameter float c          := 1.0
parameter float epsilon    := 1.0

parameter float c_pi       := 3.14159265358979324

# Coordinate system, grid 
cartesian coordinates t, x

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

# Grid functions
float v  on g1 at 0,1  {out_gf := 1}
float w  on g1 at 0,1  {out_gf := 1}

# Difference operators 
operator DP(f,t)       := (<1>f[0] -  <0>f[0]) / dt
operator MU(f,t)       := (<1>f[0] + <0>f[0]) / 2
operator D0(f,x)       := (<0>f[1] - <0>f[-1])  / (2 * dx)
operator DP(f,x)       := (<0>f[1] - <0>f[0])  / dx
operator DM(f,x)       := (<0>f[0] - <0>f[-1]) / dx

# Difference equations 
evaluate residual v { 
   [1:1]    := <1>v[0] = 0;
   [2:Nx-1] := DP(v,t) = MU(D0(w,x) + (-v + c * w) / epsilon,t);
   [Nx:Nx]  := <1>v[0] = 0;
}

evaluate residual w { 
   [1:1]    := DP(w,t) = MU(DP(v,x),t);
   [2:Nx-1] := DP(w,t) = MU(D0(v,x),t);
   [Nx:Nx]  := DP(w,t) = MU(DM(v,x),t);
}

initialize v {
   [1:Nx] := amp * exp(-((x-x0)/delta)^2);
}

initialize w {
   [1:Nx] := (v + exp(-(x-0.5)^2/0.01) * epsilon) / c;
}

looper iterative

auto initialize v, w

auto update v, w