########################################################################
# Project 2 example: Klein-Gordon equation in spherically symmetric    #
#                    Minkowski spacetime, demonstrating how            #
#                    dissipation is added                              #
########################################################################

########################################################################
# Parameters                                                           #
########################################################################

# This is how to set the memory size
system parameter int memsiz := 20000000

########################################################################
#                                                                      #
# Initial pulse parameters:                                            #
#                                                                      #
# let phi(r,t=0) = A0 * exp(-(r-r0)^2/Delta^2)                         #
#                                                                      #
# then xi(r,t=0) = d(phi)/dr                                           #
#                                                                      #
# and choose PI(r,t=0) = xi(r,t=0)                                     #
#                                                                      #
# for an approximately ingoing wave at t=0 (assuming r0>>Delta)        #
#                                                                      #
########################################################################

parameter float A0 := 1
parameter float r0 := 10
parameter float Delta := 1

########################################################################
# Kreiss-Oliger dissipation parameter ( 0 <= eps < 1)                  #
########################################################################

parameter float eps := 0.5

########################################################################
# Coordinate definition                                                # 
########################################################################
rect coordinates t,r

parameter float rmin := 0.0
parameter float rmax := 20.0

uniform rect grid g1 [1:Nr] {rmin:rmax}

########################################################################
# Use a 2-time level, Crank-Nicholson style differencing scheme        #
########################################################################
operator AVG(f,t) := (<1>f[0]+<0>f[0])/2
operator FW(f,t) := (<1>f[0])
operator D_LF(f,r) := (<0>f[1] - <0>f[-1])/(2*dr)
operator D_FW(f,r) := (-3*<0>f[0]+4*<0>f[1]-<0>f[2])/(2*dr)
operator D_BW(f,r) := (3*<0>f[0]-4*<0>f[-1]+<0>f[-2])/(2*dr)

operator D_CN(f,t) := (<1>f[0] - <0>f[0])/(dt)
operator D_CN(f,r) := AVG(D_LF(f,r),t)
operator D_CNFW(f,r) := AVG(D_FW(f,r),t)
operator D_CNBW(f,r) := AVG(D_BW(f,r),t)

########################################################################
# Kreiss-Oliger dissipation operator                                   #
########################################################################
operator DISS_KO(f,t) := -(eps/(16*dt))* (<0>f[-2]-4*<0>f[-1]+6*<0>f[0]-4*<0>f[1]+<0>f[2])
operator DISS_KO_ODD_FW1(f,t) := -(eps/(16*dt))* (-<0>f[0]-4*<0>f[-1]+6*<0>f[0]-4*<0>f[1]+<0>f[2])
operator DISS_KO_EVEN_FW1(f,t) := -(eps/(16*dt))* (<0>f[0]-4*<0>f[-1]+6*<0>f[0]-4*<0>f[1]+<0>f[2])

########################################################################
# Variables                                                            #
########################################################################
float PI    on g1 at 0,1 { out_gf := 1 }
float xi    on g1 at 0,1 { out_gf := 1 }

########################################################################
# Initial conditions                                                   #
########################################################################

initialize xi 
{
   [1:Nr] := -2*(r-r0)/Delta^2 * A0*exp(-(r-r0)^2/Delta^2);
}

initialize PI 
{
   [1:Nr] := -2*(r-r0)/Delta^2 * A0*exp(-(r-r0)^2/Delta^2);
}

########################################################################
# Residuals (equations of motion)                                      #
########################################################################

evaluate residual xi
{
   [1:1]       := <1>xi[0]=0;
   [2:2]       := D_CN(xi,t) = D_CN(PI,r) + DISS_KO_ODD_FW1(xi,t);
   [3:Nr-2]    := D_CN(xi,t) = D_CN(PI,r) + DISS_KO(xi,t);
   [Nr-1:Nr-1] := D_CN(xi,t) = D_CN(PI,r);
   [Nr:Nr]     := D_CN(xi,t) + AVG(xi/r+D_BW(xi,r),t) = 0;
}

evaluate residual PI
{
   [1:1]       := D_CNFW(PI,r) = 0;
   [2:2]       := D_CN(PI,t) = AVG(1/r^2,t)*D_CN(r^2*xi,r) + DISS_KO_EVEN_FW1(PI,t);
   [3:Nr-2]    := D_CN(PI,t) = AVG(1/r^2,t)*D_CN(r^2*xi,r) + DISS_KO(PI,t);
   [Nr-1:Nr-1] := D_CN(PI,t) = AVG(1/r^2,t)*D_CN(r^2*xi,r);
   [Nr:Nr]     := D_CN(PI,t) + AVG(PI/r+D_BW(PI,r),t) = 0;
}

########################################################################
# Define functions to initialize (note order)                          #
########################################################################

auto initialize PI,xi

########################################################################
# Define order to evolve, and solve via iteration                      #
########################################################################

auto update xi,PI

looper iterative