## import packages
using MAT,Plots,Dates,TimerOutputs,WriteVTK,DataFrames,CSV,
Statistics,DelimitedFiles,ProgressMeter

const USE_GPU=false;  # Use GPU? If this is set false, then no GPU needs to be available
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
    @init_parallel_stencil(CUDA,Float64,3);
else
    @init_parallel_stencil(Threads,Float64,3);
end
## 8th-order derivatives
@parallel function Dx_inn(in,out)
    @inn(out)=@d_xi_8(in);
    return nothing
end

@parallel function Dy_inn(in,out)
    @inn(out)=@d_yi_8(in);
    return nothing
end

@parallel function Dz_inn(in,out)
    @inn(out)=@d_zi_8(in);
    return nothing
end

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

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

@parallel_indices (ix,iz) function u_2_plus(in,out)
out[ix,4:end-2,iz]=in[ix,:,iz];
return nothing
end

@parallel_indices (ix,iz) function u_2_minus(in,out)
out[ix,3:end-3,iz]=in[ix,:,iz];
return nothing
end

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

@parallel_indices (ix,iy) function u_3_minus(in,out)
out[ix,iy,3:end-3]=in[ix,iy,:];
return nothing
end
##
nx=12;
ny=12;
nz=12;
A=@zeros(nx,ny,nz);
dtt1=@zeros(nx-5,ny,nz);
dtt2=@zeros(nx,ny-5,nz);
dtt3=@zeros(nx,ny,nz-5);
C=@zeros(nx,ny,nz);
A[5,5,5]=1;
## Dx_plus: place derivative in the later point
@parallel Dx_inn(A,dtt1);
@parallel (2:ny-1,2:nz-1) u_1_plus(dtt1,C);
## Dx_minus: place derivative in the previous point
@parallel Dx_inn(A,dtt1);
@parallel (2:ny-1,2:nz-1) u_1_minus(dtt1,C);
## Dy_plus: place derivative in the later point
@parallel Dy_inn(A,dtt2);
@parallel (2:nx-1,2:nz-1) u_2_plus(dtt2,C);
## Dy_minus: place derivative in the previous point
@parallel Dy_inn(A,dtt2);
@parallel (2:nx-1,2:nz-1) u_2_minus(dtt2,C);
## Dz_plus: place derivative in the later point
@parallel Dz_inn(A,dtt3);
@parallel (2:nx-1,2:ny-1) u_3_plus(dtt3,C);
## Dz_minus: place derivative in the previous point
@parallel Dz_inn(A,dtt3);
@parallel (2:nx-1,2:nz-1) u_3_minus(dtt3,C);
## 2nd-order derivatives
@parallel function Dx_inn(in,out)
    @inn(out)=@d_xi_1(in);
    return nothing
end

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

@parallel_indices (iy,iz) function u_1_minus(in,out)
out[1:end-2,iy,iz]=in[:,iy,iz];
return nothing
end
##
nx=12;
ny=12;
nz=12;
IND_accuracy=2;
A=@zeros(nx,ny,nz);
dtt1=@zeros(nx-IND_accuracy,ny,nz);
dtt2=@zeros(nx,ny-IND_accuracy,nz);
dtt3=@zeros(nx,ny,nz-IND_accuracy);
C=@zeros(nx,ny,nz);
A[5,5,5]=1;
## Dx_plus: place derivative in the later point
@parallel Dx_inn(A,dtt1);
@parallel (2:ny-1,2:nz-1) u_1_plus(dtt1,C);
## Dx_minus: place derivative in the previous point
@parallel Dx_inn(A,dtt1);
@parallel (2:ny-1,2:nz-1) u_1_minus(dtt1,C);
## 12th-order derivatives
@parallel function Dx_inn(in,out)
    @inn(out)=@d_xi_12(in);
    return nothing
end

@parallel function Dy_inn(in,out)
    @inn(out)=@d_yi_12(in);
    return nothing
end

@parallel function Dz_inn(in,out)
    @inn(out)=@d_zi_12(in);
    return nothing
end

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

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

@parallel_indices (ix,iz) function u_2_plus(in,out)
out[ix,6:end-4,iz]=in[ix,:,iz];
return nothing
end

@parallel_indices (ix,iz) function u_2_minus(in,out)
out[ix,5:end-5,iz]=in[ix,:,iz];
return nothing
end

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

@parallel_indices (ix,iy) function u_3_minus(in,out)
out[ix,iy,5:end-5]=in[ix,iy,:];
return nothing
end
##
nx=21;
ny=21;
nz=21;
A=@zeros(nx,ny,nz);
dtt1=@zeros(nx-9,ny,nz);
dtt2=@zeros(nx,ny-9,nz);
dtt3=@zeros(nx,ny,nz-9);
C=@zeros(nx,ny,nz);
A[11,11,11]=1;
## Dx_plus: place derivative in the later point
@parallel Dx_inn(A,dtt1);
@parallel (2:ny-1,2:nz-1) u_1_plus(dtt1,C);
## Dx_minus: place derivative in the previous point
@parallel Dx_inn(A,dtt1);
@parallel (2:ny-1,2:nz-1) u_1_minus(dtt1,C);
