# RNPL program for solution of wavemap as described by Lars Andersson, 2002/02

# Interior Crank-Nicholson scheme with O(h^2) centred spatial derivatives.
# Outgoing boundary conditions based on approximate 1/sqrt(r) fall-off

# 2007-08-08:  Implementing maintenance of constraints via simple reprojection
#    onto the unit sphere

# 2007-08-27:  Implementing equivariant initial data winding number = 1

# Set working memory size
system parameter int memsiz := 160000000

# Computing 
parameter float xmin := -1
parameter float xmax := 1
parameter float ymin := -1
parameter float ymax := 1

parameter float epsilon := -1.0

parameter float amp := 0.3
parameter float cx := 0.0
parameter float cy := 0.0
parameter float sx := 0.1
parameter float sy := 0.1

parameter float cr   := 0.5
parameter float sr   := 0.2

parameter float amp1 := 0.5
parameter float cr1  := 0.5
parameter float sr1  := 0.2
parameter float cy1  := 0.0

parameter float amp2 := 0.5
parameter float cr2  := 0.3
parameter float sr2  := 0.2
parameter float cy2  := 1.0

parameter int   consrp := 0

# Current computed maximum of source function
parameter float s_max       := 0

# parameters for parallel code ... if one, the corresponding boundary
# is a physical boundary.
parameter int xmax_pb := 1
parameter int xmin_pb := 1
parameter int ymax_pb := 1
parameter int ymin_pb := 1

rec coordinates t,x,y

uniform rec grid g1 [1:Nx][1:Ny] {xmin:xmax}{ymin:ymax} 
attribute int out_gf encodeone := [ 0 ]

# 3-vector of wave fields
float u1  on g1 at 0, 1 {out_gf := 1}
float u2  on g1 at 0, 1 {out_gf := 1}
float u3  on g1 at 0, 1 {out_gf := 1}

# Field time derivatives
float pi1 on g1 at 0, 1 {out_gf := 0}
float pi2 on g1 at 0, 1 {out_gf := 0}
float pi3 on g1 at 0, 1 {out_gf := 0}

float temp1 on g1 at 0, 1 

# Nonlinear source function
float s on g1 at 0, 1 {out_gf := 1}

# r, sqrt(r)
float r  on g1
float rh on g1

# Auxiliary functions for initial data
float phi on g1       {out_gf := 0}
float sth on g1       {out_gf := 0}
float cth on g1       {out_gf := 0}

# Grid-function for checking preservation of "constraint" 
# u1^2 + u2^2 + u3^2 = 1
float cons on g1  {out_gf := 0}

attribute int out_gf encodeone := [ 0 ]

# Forward time shift operator
operator F(f,t)    := <1>f[0][0]
# Crank-Nicholson time derivative
operator DCN(f,t)  := (<1>f[0][0] - <0>f[0][0])/dt
# Time averaging operator
operator MU(f,t)    := (<1>f[0][0] + <0>f[0][0])/2
# Kreiss-Oliger dissipation operator
operator DISS(f,t)  := -epsdis / (16*dt) * (12 *<0>f[0][0] + 
    <0>f[2][0] + <0>f[-2][0] - 4 * (<0>f[1][0] + <0>f[-1][0]) +
    <0>f[0][2] + <0>f[0][-2] - 4 * (<0>f[0][1] + <0>f[0][-1]))

# Centred first and second spatial derivatives
operator D0(f,x)   := (<0>f[1][0] - <0>f[-1][0])/(2*dx)
operator D0(f,x,x) := (<0>f[1][0] - 2*<0>f[0][0] + <0>f[-1][0])/(dx^2)
operator D0(f,y)   := (<0>f[0][1] - <0>f[0][-1])/(2*dy)
operator D0(f,y,y) := (<0>f[0][1] - 2*<0>f[0][0] + <0>f[0][-1])/(dy^2)

# Backwards/forwards second-order first spatial derivatives 
operator DB(f,x)   := ( 3*<0>f[0][0] - 4*<0>f[-1][0] + <0>f[-2][0])/(2*dx)
operator DF(f,x)   := (-3*<0>f[0][0] + 4*<0>f[1][0] - <0>f[2][0])/(2*dx)
operator DB(f,y)   := ( 3*<0>f[0][0] - 4*<0>f[0][-1] + <0>f[0][-2])/(2*dy)
operator DF(f,y)   := (-3*<0>f[0][0] + 4*<0>f[0][1] - <0>f[0][2])/(2*dy)

# Definition of source function
residual s {
          [1:Nx][1:1] := if ( ymin_pb == 1 ) then F(s,t) = 0;
        [1:Nx][Ny:Ny] := if ( ymax_pb == 1 ) then F(s,t) = 0;
          [1:1][1:Ny] := if ( xmin_pb == 1 ) then F(s,t) = 0;
        [Nx:Nx][1:Ny] := if ( xmax_pb == 1 ) then F(s,t) = 0;
     [2:Nx-1][2:Ny-1] := F(s,t) = F((-pi1^2 + D0(u1,x)^2 + D0(u1,y)^2)  +
                             (-pi2^2 + D0(u2,x)^2 + D0(u2,y)^2)  +
                             (-pi3^2 + D0(u3,x)^2 + D0(u3,y)^2),t);
}

# Equations of motion, note that the preponderance of residual 
# definitions are, as usual, for the treatment of boundary 
# conditions.
evaluate residual u1 { 
          [1:Nx][1:1] := if ( ymin_pb == 1 ) then DCN(u1,t) = MU(pi1,t);
        [1:Nx][Ny:Ny] := if ( ymax_pb == 1 ) then DCN(u1,t) = MU(pi1,t);
          [1:1][1:Ny] := if ( xmin_pb == 1 ) then DCN(u1,t) = MU(pi1,t);
        [Nx:Nx][1:Ny] := if ( xmax_pb == 1 ) then DCN(u1,t) = MU(pi1,t);

     [2:Nx-1][2:Ny-1] := DCN(u1,t) = MU(pi1,t);
}

evaluate residual pi1 { 
        [2:Nx-1][1:1] := if( ymin_pb == 1 ) then DCN(rh*pi1,t) + 
                           MU((x/r)*D0(rh*pi1,x) + (y/r)*DF(rh*pi1,y),t) = 0;
      [2:Nx-1][Ny:Ny] := if( ymax_pb == 1 ) then DCN(rh*pi1,t) +
                           MU((x/r)*D0(rh*pi1,x) + (y/r)*DB(rh*pi1,y),t) = 0;
        [1:1][2:Ny-1] := if( xmin_pb == 1) then DCN(rh*pi1,t) +
                           MU((x/r)*DF(rh*pi1,x) + (y/r)*D0(rh*pi1,y),t) = 0;
      [Nx:Nx][2:Ny-1] := if( xmax_pb == 1) then DCN(rh*pi1,t) +
                           MU((x/r)*DB(rh*pi1,x) + (y/r)*D0(rh*pi1,y),t) = 0;
           [1:1][1:1] := if( xmin_pb == 1 && ymin_pb == 1 ) then DCN(rh*pi1,t) +
                           MU((x/r)*DF(rh*pi1,x) + (y/r)*DF(rh*pi1,y),t) = 0;
         [Nx:Nx][1:1] := if( xmax_pb == 1 && ymin_pb == 1 ) then DCN(rh*pi1,t) +
                           MU((x/r)*DB(rh*pi1,x) + (y/r)*DF(rh*pi1,y),t) = 0;
         [1:1][Ny:Ny] := if( xmin_pb == 1 && ymax_pb == 1 ) then DCN(rh*pi1,t) +
                           MU((x/r)*DF(rh*pi1,x) + (y/r)*DB(rh*pi1,y),t) = 0;
       [Ny:Ny][Ny:Ny] := if( xmax_pb == 1 && ymax_pb == 1 ) then DCN(rh*pi1,t) +
                           MU((x/r)*DB(rh*pi1,x) + (y/r)*DB(rh*pi1,y),t) = 0;

     [2:Nx-1][2:Ny-1] := DCN(pi1,t) = MU(D0(u1,x,x) + D0(u1,y,y) -
                                       epsilon*u1*s,t);
}

evaluate residual u2 { 
          [1:Nx][1:1] := if ( ymin_pb == 1 ) then DCN(u2,t) = MU(pi2,t);
        [1:Nx][Ny:Ny] := if ( ymax_pb == 1 ) then DCN(u2,t) = MU(pi2,t);
          [1:1][1:Ny] := if ( xmin_pb == 1 ) then DCN(u2,t) = MU(pi2,t);
        [Nx:Nx][1:Ny] := if ( xmax_pb == 1 ) then DCN(u2,t) = MU(pi2,t);

     [2:Nx-1][2:Ny-1] := DCN(u2,t) = MU(pi2,t);
}

evaluate residual pi2 { 
        [2:Nx-1][1:1] := if( ymin_pb == 1 ) then DCN(rh*pi2,t) + 
                           MU((x/r)*D0(rh*pi2,x) + (y/r)*DF(rh*pi2,y),t) = 0;
      [2:Nx-1][Ny:Ny] := if( ymax_pb == 1 ) then DCN(rh*pi2,t) +
                           MU((x/r)*D0(rh*pi2,x) + (y/r)*DB(rh*pi2,y),t) = 0;
        [1:1][2:Ny-1] := if( xmin_pb == 1) then DCN(rh*pi2,t) +
                           MU((x/r)*DF(rh*pi2,x) + (y/r)*D0(rh*pi2,y),t) = 0;
      [Nx:Nx][2:Ny-1] := if( xmax_pb == 1) then DCN(rh*pi2,t) +
                           MU((x/r)*DB(rh*pi2,x) + (y/r)*D0(rh*pi2,y),t) = 0;
           [1:1][1:1] := if( xmin_pb == 1 && ymin_pb == 1 ) then DCN(rh*pi2,t) +
                           MU((x/r)*DF(rh*pi2,x) + (y/r)*DF(rh*pi2,y),t) = 0;
         [Nx:Nx][1:1] := if( xmax_pb == 1 && ymin_pb == 1 ) then DCN(rh*pi2,t) +
                           MU((x/r)*DB(rh*pi2,x) + (y/r)*DF(rh*pi2,y),t) = 0;
         [1:1][Ny:Ny] := if( xmin_pb == 1 && ymax_pb == 1 ) then DCN(rh*pi2,t) +
                           MU((x/r)*DF(rh*pi2,x) + (y/r)*DB(rh*pi2,y),t) = 0;
       [Ny:Ny][Ny:Ny] := if( xmax_pb == 1 && ymax_pb == 1 ) then DCN(rh*pi2,t) +
                           MU((x/r)*DB(rh*pi2,x) + (y/r)*DB(rh*pi2,y),t) = 0;

     [2:Nx-1][2:Ny-1] := DCN(pi2,t) = MU(D0(u2,x,x) + D0(u2,y,y) -
                                       epsilon*u2*s,t);
}

evaluate residual u3 { 
          [1:Nx][1:1] := if ( ymin_pb == 1 ) then DCN(u3,t) = MU(pi3,t);
        [1:Nx][Ny:Ny] := if ( ymax_pb == 1 ) then DCN(u3,t) = MU(pi3,t);
          [1:1][1:Ny] := if ( xmin_pb == 1 ) then DCN(u3,t) = MU(pi3,t);
        [Nx:Nx][1:Ny] := if ( xmax_pb == 1 ) then DCN(u3,t) = MU(pi3,t);

     [2:Nx-1][2:Ny-1] := DCN(u3,t) = MU(pi3,t);
}

evaluate residual pi3 { 
        [2:Nx-1][1:1] := if( ymin_pb == 1 ) then DCN(rh*pi3,t) + 
                           MU((x/r)*D0(rh*pi3,x) + (y/r)*DF(rh*pi3,y),t) = 0;
      [2:Nx-1][Ny:Ny] := if( ymax_pb == 1 ) then DCN(rh*pi3,t) +
                           MU((x/r)*D0(rh*pi3,x) + (y/r)*DB(rh*pi3,y),t) = 0;
        [1:1][2:Ny-1] := if( xmin_pb == 1) then DCN(rh*pi3,t) +
                           MU((x/r)*DF(rh*pi3,x) + (y/r)*D0(rh*pi3,y),t) = 0;
      [Nx:Nx][2:Ny-1] := if( xmax_pb == 1) then DCN(rh*pi3,t) +
                           MU((x/r)*DB(rh*pi3,x) + (y/r)*D0(rh*pi3,y),t) = 0;
           [1:1][1:1] := if( xmin_pb == 1 && ymin_pb == 1 ) then DCN(rh*pi3,t) +
                           MU((x/r)*DF(rh*pi3,x) + (y/r)*DF(rh*pi3,y),t) = 0;
         [Nx:Nx][1:1] := if( xmax_pb == 1 && ymin_pb == 1 ) then DCN(rh*pi3,t) +
                           MU((x/r)*DB(rh*pi3,x) + (y/r)*DF(rh*pi3,y),t) = 0;
         [1:1][Ny:Ny] := if( xmin_pb == 1 && ymax_pb == 1 ) then DCN(rh*pi3,t) +
                           MU((x/r)*DF(rh*pi3,x) + (y/r)*DB(rh*pi3,y),t) = 0;
       [Ny:Ny][Ny:Ny] := if( xmax_pb == 1 && ymax_pb == 1 ) then DCN(rh*pi3,t) +
                           MU((x/r)*DB(rh*pi3,x) + (y/r)*DB(rh*pi3,y),t) = 0;

     [2:Nx-1][2:Ny-1] := DCN(pi3,t) = MU(D0(u3,x,x) + D0(u3,y,y) -
                                       epsilon*u3*s,t);
}

residual cons {
        [1:Nx][1:Ny]  := cons = u1^2 + u2^2 + u3^2 - 1;
}


#############################################################################
# Initial data for time symmetric, distorted, gaussian rings in u1, u2
#############################################################################
initialize r {
   [1:Nx][1:Ny]   := sqrt(x^2 + y^2);
}

initialize rh {
   [1:Nx][1:Ny]   := sqrt(sqrt(x^2 + y^2));
}

initialize sth {
   [1:Nx][1:Ny]   := if x == 0 && y == 0 then 0 else y / r;
}

initialize cth {
   [1:Nx][1:Ny]   := if x == 0 && y == 0 then 0 else x / r;
}

initialize phi {
   [1:Nx][1:Ny] := amp * exp(-((r - cr) / sr)^2);
}

initialize u1 { 
   [1:Nx][1:Ny] := sin(phi) * cth; 
}
initialize u2 { 
   [1:Nx][1:Ny] := sin(phi) * sth;
}
initialize u3 { 
   [1:Nx][1:Ny] := cos(phi);
}

initialize pi1 { 
   [1:Nx][1:Ny]   := 0;
}

initialize pi2 { 
   [1:Nx][1:Ny]   := 0;
}
initialize pi3 { 
   [1:Nx][1:Ny]   := 0;
}
#############################################################################

initialize s {
         [1:Nx][1:1]  := 0;
        [1:Nx][Ny:Ny] := 0;
          [1:1][1:Ny] := 0;
        [Nx:Nx][1:Ny] := 0;
     [2:Nx-1][2:Ny-1] := (-pi1^2 + D0(u1,x)^2 + D0(u1,y)^2)  +
                         (-pi2^2 + D0(u2,x)^2 + D0(u2,y)^2)  +
                         (-pi3^2 + D0(u3,x)^2 + D0(u3,y)^2);
}

looper iterative

auto update u1, u2, u3, pi1, pi2, pi3, s, cons

auto initialize r, rh, cth, sth, phi
auto initialize u1, u2, u3, pi1, pi2, pi3, s

cons_rp.inc cons_rp UPDATES temp1 HEADER u1, u2, u3, s, consrp, t, 
	s_max