############################################################
# Heat equation
# 
#    u_t - u_xx = 0
#
# on 0 <= x <= 1 with 
#
#    u(0) = u(1) = 0
#
# Crank-Nicholson differencing, multigrid solution
############################################################

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

parameter float amp        := 1.0

parameter float xmin       := 0.0
parameter float xmax       := 1.0
parameter float xc         := 0.5
parameter float xwid       := 0.05

parameter float ymin       := 0.0
parameter float ymax       := 1.0
parameter float yc         := 0.5
parameter float ywid       := 0.05

parameter float mpi        := 3.1415926535897932
parameter float camp       := 0.0

parameter float mgeps      := 1.0e-9
parameter int   mgmaxcyc   := 50

rect coordinates t, x, y

uniform rect grid g1 [1:Nx][1:Ny] {xmin:xmax}{ymin:ymax}

float u    on g1 at 0,1 {out_gf := 1}
float rhs  on g1 
float ures on g1 at 0,1 {out_gf := 1}

# Crank Nicholson time derivative (first forward difference)
operator D(f,t) := (<1>f[0][0] - <0>f[0][0]) / dt

# Forward averaging operator
operator MU(f,t) := (<1>f[0][0] + <0>f[0][0]) / 2

# O(h^2) centred and backward differences 
operator D0(f,x,x) := (<0>f[1][0] - 2*<0>f[0][0] + <0>f[-1][0]) / (dx^2)
operator D0(f,y,y) := (<0>f[0][1] - 2*<0>f[0][0] + <0>f[0][-1]) / (dy^2)

# Difference equations
residual ures { 
            [2:Nx-1][2:Ny-1] := ures = D(u,t)- MU(D0(u,x,x),t) - MU(D0(u,y,y),t);
               [1:1][1:Ny]   := ures = 0;
             [Nx:Nx][1:Ny]   := ures = 0;
              [1:Nx][1:1]    := ures = 0;
              [1:Nx][Ny:Ny]  := ures = 0;
}

# Initialize to gaussian pulse
initialize u {[1:Nx][1:Ny] := amp*exp(-((x-xc)^2/xwid^2 + (y-yc)^2/ywid^2))}
initialize ures {[1:Nx][1:Ny] := 0}

looper standard

update0.inc update0 updates u header u[unp1,un], rhs, x, y, t, dx, dy, dt, mgeps, mgmaxcyc
auto update ures