############################################################
# Equivariant 2+1D wave map
# 
#    ph_tt = ((r ph_r)_r) / r - sin(2 ph) / 2 r^2
#
# solved in second order form 
#
#    ph_t  = pi
#    pi_t  = ((r ph_r)_r) / r - sin(2 ph) / 2 r^2
#
# for integer p >= 1
# 
# Uses Crank-Nicholson with (implicit) Kreiss-Oliger 
# dissipation and (implicit) ingoing/outgoing radiation 
# conditions (O(h^2) forwards and backwards differences) at 
# the boundaries.  
# 
# Initial data: 
# 
#    Time symmetric gaussian
#
############################################################

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

parameter float rmin       :=  0.0
parameter float rmax       := 32.0

# K/O dissipation parameter
parameter float epsdis     := 0.5

# Coupling and power parameters 
parameter float kappa      := 1.0
parameter int   p          := 1

parameter float amp        := 1.0e-1
parameter float r0         := 0.0
parameter float rwid       := 1.0
parameter float idsignum   := 1.0

# For interface to adaptive framework
parameter int   swlb       := 1
parameter int   swrb       := 1

spherical coordinates t, r

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

float ph   on g1 at 0,1 {out_gf := 1}
float pi   on g1 at 0,1 {out_gf := 1}
float rh   on g1

# Time-translation operator
operator T(f,t) := <1>f[0]

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

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

operator MUM(f,r) := (<0>f[0] + <0>f[-1]) / 2

# Dissipation operator
operator DISS(f,t) := -epsdis / (16 * dt) *
        (6 *<0>f[0] + <0>f[2] + <0>f[-2] - 4 * (<0>f[1] + <0>f[-1]))

# O(h^2) centred, forward and backward first derivatives
operator D0(f,r)   := (<0>f[1] - <0>f[-1]) / (2 * dr)
operator D0(f,r,r) := (<0>f[1] - 2*<0>f[0] + <0>f[-1]) / dr^2

# Compact, origin-regularized differencing of radial part of Laplacian
operator D0R(f,r,r) := ((r + 0.5*dr) * (<0>f[1] - <0>f[0]) - 
                        (r - 0.5*dr) * (<0>f[0] - <0>f[-1])) /
                       (r * dr^2)

operator DM(f,r) := (<0>f[0] - <0>f[-1]) / dr
operator DB(f,r) := ( 3*<0>f[0] - 4*<0>f[-1] + <0>f[-2]) / (2 * dr)
operator DF(f,r) := (-3*<0>f[0] + 4*<0>f[1]  - <0>f[2])  / (2 * dr)

# Difference equations

evaluate residual ph { 
               [1:1] := if ( swlb == 1 ) then T(ph,t) = 0;
               [2:2] := D(ph,t)  =  MU(pi,t); 
            [3:Nr-2] := D(ph,t)  =  MU(pi + DISS(ph,t),t);
         [Nr-1:Nr-1] := D(ph,t)  =  MU(pi,t); 
             [Nr:Nr] := if ( swrb == 1 ) then D(rh* ph,t)  =  MU(-DB(rh * ph,r),t);
}

evaluate residual pi { 
               [1:1] := if ( swlb == 1 ) then T(pi,t) = 0;
               [2:2] := D(pi,t) =  MU(D0R(ph,r,r) - kappa * sin(2*ph)/(2*r^2),t);
            [3:Nr-2] := D(pi,t) =  MU(D0R(ph,r,r) - kappa * sin(2*ph)/(2*r^2) + DISS(pi,t),t);
         [Nr-1:Nr-1] := D(pi,t) =  MU(D0R(ph,r,r) - kappa * sin(2*ph)/(2*r^2),t);
             [Nr:Nr] := if ( swrb == 1 ) then D(rh * pi,t) =  MU(-DB(rh * pi,r),t);
}

initialize rh {
	[1:Nr] :=  sqrt(r);
}

initialize ph {
	[1:1]  := if ( swlb == 1 ) then 0 else amp * exp(-((r - r0) / rwid)^2 );
	[2:Nr] :=  amp * exp(-((r - r0) / rwid)^2 );
}

initialize pi {
	[1:1]    := if ( swlb == 1 ) then 0 else idsignum * D0(rh * ph,r) / rh;
	[2:Nr-1] := idsignum * D0(rh * ph,r) / rh;
	[Nr:Nr]  := idsignum * DB(rh * ph,r) / rh;
}

looper iterative

auto initialize rh
auto initialize ph, pi
auto update ph, pi