## This is to test the first derivative of sin
using MAT,Plots,LinearAlgebra
const USE_GPU=false;  # Use GPU? If this is set false, then no GPU needs to be available
using ParallelStencil
using ParallelStencil.FiniteDifferences1D
@static if USE_GPU
    @init_parallel_stencil(CUDA,Float64,3);
else
    @init_parallel_stencil(Threads,Float64,3);
end
dx=.001;
x=0:dx:6*pi;
nx=length(x);
y=cos.(x);
## analytical expression of d^2y/dx^2
d2y_dx2=-cos.(x);
## d^2y/dx^2 2nd-order accuracy
@parallel function Dx_inn_2nd_order_accuracy(in,out,dx)
    @inn(out)=@d_1(in)/dx;
    return nothing
end

@parallel_indices (iy) function u_1_plus_2nd_order_accuracy(in,out)
out[2:end-1]=in[:];
return nothing
end

@parallel_indices (iy) function u_1_minus_2nd_order_accuracy(in,out)
out[1:end-2]=in[:];
return nothing
end

dtt=@zeros(nx-2,);
tt=@zeros(nx,);
d2y_dx2_2nd_order_accuracy=@zeros(nx,);

@parallel Dx_inn_2nd_order_accuracy(y,dtt,dx);
@parallel (1:1) u_1_minus_2nd_order_accuracy(dtt,tt);
@parallel Dx_inn_2nd_order_accuracy(tt,dtt,dx);
@parallel (1:1) u_1_plus_2nd_order_accuracy(dtt,d2y_dx2_2nd_order_accuracy);
## d^2y/dx^2 8th-order accuracy
@parallel function Dx_inn_8th_order_accuracy(in,out,dx)
    @inn(out)=@d_8(in)/dx;
    return nothing
end

@parallel_indices (iy) function u_1_plus_8th_order_accuracy(in,out)
out[4:end-2]=in[:];
return nothing
end

@parallel_indices (iy) function u_1_minus_8th_order_accuracy(in,out)
out[3:end-3]=in[:];
return nothing
end

dtt=@zeros(nx-5,);
tt=@zeros(nx,);
d2y_dx2_8th_order_accuracy=@zeros(nx,);

@parallel Dx_inn_8th_order_accuracy(y,dtt,dx);
@parallel (1:1) u_1_minus_8th_order_accuracy(dtt,tt);
@parallel Dx_inn_8th_order_accuracy(tt,dtt,dx);
@parallel (1:1) u_1_plus_8th_order_accuracy(dtt,d2y_dx2_8th_order_accuracy);
## d^2y/dx^2 12th-order accuracy
@parallel function Dx_inn_12th_order_accuracy(in,out,dx)
    @inn(out)=@d_12(in)/dx;
    return nothing
end

@parallel_indices (iy) function u_1_plus_12th_order_accuracy(in,out)
out[6:end-4]=in[:];
return nothing
end

@parallel_indices (iy) function u_1_minus_12th_order_accuracy(in,out)
out[5:end-5]=in[:];
return nothing
end

dtt=@zeros(nx-9,);
tt=@zeros(nx,);
d2y_dx2_12th_order_accuracy=@zeros(nx,);

@parallel Dx_inn_12th_order_accuracy(y,dtt,dx);
@parallel (1:1) u_1_minus_12th_order_accuracy(dtt,tt);
@parallel Dx_inn_12th_order_accuracy(tt,dtt,dx);
@parallel (1:1) u_1_plus_12th_order_accuracy(dtt,d2y_dx2_12th_order_accuracy);
## plot derivatives
IND=30:nx-30;

plot(x[IND],d2y_dx2[IND],label="analytical",color="red",title="d^2y/dx^2",lw=3);
plot!(x[IND],d2y_dx2_2nd_order_accuracy[IND],label="2nd-order",color="blue",lw=3)
plot!(x[IND],d2y_dx2_8th_order_accuracy[IND],label="8th-order",color="green",lw=3)
plot!(x[IND],d2y_dx2_12th_order_accuracy[IND],label="12th-order",color="yellow",lw=3)
## plot error
plot(x[IND],d2y_dx2_2nd_order_accuracy[IND]-d2y_dx2[IND],label="2nd-order",color="blue",lw=3);
plot!(x[IND],d2y_dx2_8th_order_accuracy[IND]-d2y_dx2[IND],label="8th-order",color="green",lw=3);
plot!(x[IND],d2y_dx2_12th_order_accuracy[IND]-d2y_dx2[IND],label="12th-order",color="yellow",lw=3)
## norm error
e1=norm(d2y_dx2_2nd_order_accuracy[IND]-d2y_dx2[IND],2)
e2=norm(d2y_dx2_8th_order_accuracy[IND]-d2y_dx2[IND],2)
e3=norm(d2y_dx2_12th_order_accuracy[IND]-d2y_dx2[IND],2)
