!Traveling taylor-green with forcing
!                                      By Taku & Pietro
! 1. PROBLEM
! 1-1. Basic equations are as follows:
!
!                          u_x +v_y=0                                   (1) 
!            u_t + u u_x + v u_y + P_x = nu (u_xx + u_yy) +F,           (2)
!            v_t + u v_x + v v_y + P_y = nu (v_xx + v_yy) +G,           (3)
!
!
!             F=\sin(x-u_0 t)\cos(y-v_0 t)[2\nu \cos t-\sin t],         (4)
!             G=-\cos(x-u_0 t)\sin(y-v_0 t)[2\nu \cos t-\sin t]         (5)
!
! The domain is (0<x<2 Pi) x (0<y<2 Pi) and the periodic boundary comdition is considered.
! The initial condition corresponds to the following exact solution.
!   
!   u=u_0+sin(x-u_0 t)cos(y-v_0 t)cos t,
!   v=v_0-cos(x-u_0 t)sin(y-v_0 t)cos t,
!   P=(1/4)[cos(2(x-u_0 t))+cos(2(y-v_0 t))]cos^2 t.
!
!
! 2. New Artificial compressibility Method
!
!      We replace the solenoidalcondition (1) by
!
!                 beta eps^2 (P_t + gamma P)+ u_x +v_y=0,                (6) 
! where eps is the lattice spacing and beta and gamma are the positive constants 
! of the order of unity. The term gamma P acts as a damping term for the acoustic mode. 
! The momentum equations (2) and (3) are left as they are (the bulk viscosity is not introduced).
!
!
!  The spatial derivatives are approximated by the conventional 3 point finite difference formulas,
!  the explicit definition of which are shown in the last part of the code. 

! We employ the semi implicit second order accurate RK method for the time integration:
!
!
!   u(t=(n+1/2)dt)=u(t=n dt)+dt/2 *(du/dt)(t=n dt)
!   v(t=(n+1/2)dt)=v(t=n dt)+dt/2 *(dv/dt)(t=n dt)
!   P(t=(n+1/2)dt)=(P(t=n dt)-dt/2 div(u,v)(t=(n+1/2))/beta/eps^2)/(1+dt/2 gamma)
!
!   u(t=(n+1)dt)=u(t=n dt)+dt *(du/dt)(t=(n+1/2))
!   v(t=(n+1)dt)=v(t=n dt)+dt *(dv/dt)(t=(n+1/2))
!   P(t=(n+1)dt)=P(t=n dt)+dt *(-div(u,v)(t=n+1/2)dt)/beta/eps^2+checkerboard *medicine)
!
!   where checkerboard=eps*C (C is a positive constatnt of the order of unity) and
!  medicine=(Laplacian P+2 v_x u_y-u_x v_y)(t=(n+1/2)dt).
! The term checkerboard * medicine becomes O(eps^3) in the diffusive mode
! and is employed for the prevention and suppression of the checkerboard instability.
!
!  In the above formulas, the spacial derivatives except div(u,v) are approximated by
!  three point finite difference operators del_x, del_y, del_xx, and del_yy
! (see the last part of the code for the explicit definitions).
!
!  div(u,v) are u_x+v_y and approximated in the following way.
!
!  del_x u= u_x +(eps^2/6) u_xxx +O(eps^4)
!
!  In the diffusive mode, u_x+v_y=O(eps^2). Then,
!
! u_xxx =-v_xxy +O(eps^2)
!       =-del_xx del_y v + O(eps^2).
!
!   v_y is computed in the same way.
!
!
!=========================================================================================================
!******************************************************* NUMBER OF MESH POINTS ***************************
module size
integer,parameter::nx=256,ny=nx
end module
!*********************************************************************************************************

module fdcommon
real(8)::eps,beta,dt,hdt,gamma,checkerboard,med_mom
real(8)::veps,hveps,vveps,hvveps,qvveps,veps12,vveps12,hvvveps,vvvveps
real(8)::vbvveps
end module

module uvp
use size
real(8),dimension(-2:nx+1,-2:ny+1):: u,ua,v,va,P,Pa
real(8),dimension(0:nx,0:ny)::div_u,medicine
end module

module viscosity
real(8)::nu
end module


module xy
use size
real(8),dimension(0:nx)::x
real(8),dimension(0:ny)::y
end module
!
module other
use size
real(8)::ssx(0:nx),ssy(0:ny)
real(8)::pi
end module

module force
use size
real(8),dimension(0:nx,0:ny)::fx,fy
end module
!main program
!
use size
use uvp
use xy
use fdcommon
use viscosity
!use numerov_data
use other


use force

implicit none
real(8)::time,tmax
real(8)::mach,sound_speed
integer(4)::itmax,itime,isamp,i,j
real(8)::xx,yy
real(8)::uerror,verror,perror
real(8)::uexact,vexact,pexact
real(8)::uzero,vzero,t,sx,cx,sy,cy
character(1),dimension(2)::cbeta
integer::ijob
!
!************************ uzero = u_0 and vzer = v_0
uzero=0.3d0
vzero=0.6d0

pi=datan(1d0)*4d0

!****************************** nu is the viscosity (1/Re) ********************
nu=0.001d0
!****************************** eps =2 Pi/nx is the lattice spacing ***********
eps=2d0*pi/dble(nx)

call miscellaneous 


gamma=1d0
checkerboard=eps*1d0

tmax=60d0

dt=1d0*eps**2
itmax=idnint(tmax/dt)
dt=tmax/dble(itmax)
hdt=dt/2d0

ijob=1
isamp=2000
if(ijob.eq.1) then
beta=2d0
else
beta=4d0
end if



call set_fdcommon

   do i=0,nx
    x(i)=eps*dble(i)
   enddo
!
   do j=0,ny
    y(j)=eps*dble(j)
   enddo


cbeta(1)='2'
cbeta(2)='4'

if(nx.eq.32) then 
     open(1,file='32-'//cbeta(ijob)//'.dat')
else if (nx.eq.64) then
     open(1,file='64-'//cbeta(ijob)//'.dat')
else if(nx.eq.128) then
     open(1,file='128-'//cbeta(ijob)//'.dat')
else if(nx.eq.256) then
     open(1,file='256-'//cbeta(ijob)//'.dat')
	 open(2,file='256-'//cbeta(ijob)//'-history.dat')
else if(nx.eq.512) then
     open(1,file='512-'//cbeta(ijob)//'.dat')
else 
     open(1,file='1024-'//cbeta(ijob)//'.dat')
    end if

write(*,*) itmax,tmax
write(*,*) beta,beta*eps**2
write(*,*) sound_speed,dt,eps
write(*,*) checkerboard

do j=0,ny
do i=0,nx
xx=x(i)
yy=y(j)
u(i,j)=uzero+dcos(yy)*dsin(xx)
v(i,j)=vzero-dcos(xx)*dsin(yy)
P(i,j)=0.25d0*(dcos(2d0*xx)+dcos(2d0*yy))
enddo
enddo

time=0d0


do itime=1,itmax
!
t=time
do j=0,ny
do i=0,nx
sx=dsin(x(i)-uzero*t);cx=dcos(x(i)-uzero*t)
sy=dsin(y(j)-vzero*t);cy=dcos(y(j)-vzero*t)
fx(i,j)=sx*cy*(2d0*nu*dcos(t)-dsin(t))
fy(i,j)=-cx*sy*(2d0*nu*dcos(t)-dsin(t))
enddo
enddo

call periodic_u(u,v,nx,ny)
call periodic_p(p,nx,ny)

call inneruv1
call periodic_u(ua,va,nx,ny)
call divergence(ua,va,div_u)
call time_integrate_p_first
call periodic_p(pa,nx,ny)
t=time+hdt
do j=0,ny
do i=0,nx
sx=dsin(x(i)-uzero*t);cx=dcos(x(i)-uzero*t)
sy=dsin(y(j)-vzero*t);cy=dcos(y(j)-vzero*t)
fx(i,j)=sx*cy*(2d0*nu*dcos(t)-dsin(t))
fy(i,j)=-cx*sy*(2d0*nu*dcos(t)-dsin(t))
enddo
enddo

call inneruv2


call checkerboard_cure(ua,va,pa,medicine)
call time_integrate_p_second
time=dt*dble(itime);t=time
uerror=0d0
verror=0d0
perror=0d0
do j=0,ny
yy=y(j)-vzero*t
sy=dsin(yy);cy=dcos(yy)
do i=0,nx
xx=x(i)-uzero*t
sx=dsin(xx);cx=dcos(xx)
uexact=uzero+sx*cy*dcos(t)
vexact=vzero-cx*sy*dcos(t)
pexact=0.25d0*(dcos(2d0*xx)+dcos(2d0*yy))*(dcos(t))**2
!call p_shift
uerror=uerror+ssx(i)*ssy(j)*dabs(u(i,j)-uexact)
verror=verror+ssx(i)*ssy(j)*dabs(v(i,j)-vexact)
perror=perror+ssx(i)*ssy(j)*dabs(p(i,j)-pexact)
enddo
enddo
uerror=uerror/(4d0*pi**2)
verror=verror/(4d0*pi**2)
perror=perror/(4d0*pi**2)

write(2,*) time,uerror,verror,perror
write(*,*) itime,time
enddo
call p_shift
call save_data
end
!==========================================
subroutine save_data
use uvp
use xy
use force

write(1,*)'VARIABLES="x","y","u","v","P"'
write(1,*)'ZONE I=',nx+1,'J=',ny+1,'DATAPACKING=POINT'


do j=0,ny
do i=0,nx
write(1,*) x(i),y(j),u(i,j),v(i,j),p(i,j)
enddo
enddo
close(1)
return
end

!==========================================

!==========================================
subroutine inneruv1
use uvp
use xy
use fdcommon
use viscosity
use force

implicit real*8 (D)
integer::i,j
real(8)::ut,vt
real(8)::ux,uy,uxx,uyy,vx,vy,vxx,vyy,px,py
real(8):: uxy,vxy
real(8)::lap_u,lap_v


!
do j=0,ny-1
do i=0,nx-1
ux=del_x(u,i,j,nx,ny)
uxx=del_xx(u,i,j,nx,ny)
vx=del_x(v,i,j,nx,ny)
vxx=del_xx(v,i,j,nx,ny)
px=del_x(p,i,j,nx,ny)

uy=del_y(u,i,j,nx,ny)
uyy=del_yy(u,i,j,nx,ny)
vy=del_y(v,i,j,nx,ny)
vyy=del_yy(v,i,j,nx,ny)
py=del_y(p,i,j,nx,ny)

lap_u=uxx+uyy
lap_v=vxx+vyy


ut=-u(i,j)*ux-v(i,j)*uy-px+nu*lap_u+fx(i,j)
vt=-u(i,j)*vx-v(i,j)*vy-py+nu*lap_v+fy(i,j)
ua(i,j)=u(i,j)+hdt*ut
va(i,j)=v(i,j)+hdt*vt
enddo
enddo
return
end
!------------------------------------------------------------------------

!==========================================
subroutine inneruv2
use uvp
use xy
use fdcommon
use viscosity
use force

implicit real*8(D)
integer::i,j
real(8)::ut,vt
real(8)::ux,uy,uxx,uyy,vx,vy,vxx,vyy,px,py
real(8):: uxy,vxy
real(8)::lap_u,lap_v

do j=0,ny-1
do i=0,nx-1

ux=del_x(ua,i,j,nx,ny)
uxx=del_xx(ua,i,j,nx,ny)
vx=del_x(va,i,j,nx,ny)
vxx=del_xx(va,i,j,nx,ny)
px=del_x(pa,i,j,nx,ny)

uy=del_y(ua,i,j,nx,ny)
uyy=del_yy(ua,i,j,nx,ny)
vy=del_y(va,i,j,nx,ny)
vyy=del_yy(va,i,j,nx,ny)
py=del_y(pa,i,j,nx,ny)

lap_u=uxx+uyy
lap_v=vxx+vyy

ut=-ua(i,j)*ux-va(i,j)*uy-px+nu*lap_u+fx(i,j)
vt=-ua(i,j)*vx-va(i,j)*vy-py+nu*lap_v+fy(i,j)
u(i,j)=u(i,j)+dt*ut
v(i,j)=v(i,j)+dt*vt
enddo
enddo
return
end
!

!-----------------------------------------------------------------------------------
subroutine time_integrate_p_first
use uvp
use fdcommon

real(8)::pt
!
do j=0,ny-1
do i=0,nx-1
pt=-div_u(i,j)*vbvveps
pa(i,j)=(p(i,j)+pt*hdt)/(1d0+hdt*gamma)
enddo
enddo
return
end
!----------------------------------------------
subroutine time_integrate_p_second
use uvp
use fdcommon
real(8)::pt
!
do j=0,ny-1
do i=0,nx-1
pt=-div_u(i,j)*vbvveps-gamma*pa(i,j)+medicine(i,j)*checkerboard
p(i,j)=p(i,j)+pt*dt
enddo
enddo
!

return
end
!

!--------------------------------------------------------------
subroutine divergence(u,v,div_u)
use size
use fdcommon
implicit none
real(8),intent(in),dimension(-2:nx+1,-2:ny+1)::u,v
real(8),intent(out),dimension(0:nx,0:ny)::div_u
real(8)::del_x,del_y,del_xxy,del_xyy
integer::i,j

do j=0,ny-1
do i=0,nx-1
div_u(i,j)=del_x(u,i,j,nx,ny)+del_y(v,i,j,nx,ny)+&
eps**2/6d0*(del_xxy(v,i,j,nx,ny)+del_xyy(u,i,j,nx,ny))
enddo
enddo

return
end
!---------------------------------------------------------------
subroutine checkerboard_cure(u,v,p,medicine)
use size
implicit real(8)(B,D,H)
real(8),intent(in),dimension(-2:nx+1,-2:ny+1)::u,v,p
real(8),dimension(0:nx,0:ny)::medicine
integer(4)::imedicine
real(8)::ux,uy,vx,vy,dpxx,dpyy,dpss,dptt

do j=0,ny-1
do i=0,nx-1


ux=del_x(u,i,j,nx,ny);vx=del_x(v,i,j,nx,ny)
uy=del_y(u,i,j,nx,ny);vy=del_y(v,i,j,nx,ny)
dpxx=del_xx(p,i,j,nx,ny)
dpyy=del_yy(p,i,j,nx,ny)

medicine(i,j)=dpxx+dpyy+2d0*(vx*uy-ux*vy)
enddo
enddo
!
return
end
!
!****************************************************************
subroutine periodic_u(u,v,nx,ny)
real(8),dimension(-2:nx+1,-2:nx+1)::u,v
do i=0,nx-1
u(i,ny)=u(i,0)
v(i,ny)=v(i,0)

u(i,ny+1)=u(i,1)
v(i,ny+1)=v(i,1)

u(i,-1)=u(i,ny-1)
v(i,-1)=v(i,ny-1)

u(i,-2)=u(i,ny-2)
v(i,-2)=v(i,ny-2)

enddo

do j=0,ny-1
u(nx,j)=u(0,j)
v(nx,j)=v(0,j)

u(nx+1,j)=u(1,j)
v(nx+1,j)=v(1,j)

u(-1,j)=u(nx-1,j)
v(-1,j)=v(nx-1,j)

u(-2,j)=u(nx-2,j)
v(-2,j)=v(nx-2,j)

enddo


u(nx,ny)=u(0,0)
v(nx,ny)=v(0,0)

u(-1,-1)=u(nx-1,ny-1)
v(-1,-1)=v(nx-1,ny-1)

u(nx,-1)=u(0,ny-1)
v(nx,-1)=v(0,ny-1)

u(-1,ny)=u(nx-1,0)
v(-1,ny)=v(nx-1,0)

return
end

!****************************************************************
subroutine periodic_p(p,nx,ny)
real(8),dimension(-2:nx+1,-2:nx+1)::p
do i=0,nx-1
p(i,ny)=p(i,0)

p(i,ny+1)=p(i,1)

p(i,-1)=p(i,ny-1)

p(i,-2)=p(i,ny-2)

enddo

do j=0,ny-1
p(nx,j)=p(0,j)

p(nx+1,j)=p(1,j)

p(-1,j)=p(nx-1,j)

p(-2,j)=p(nx-2,j)

enddo


p(nx,ny)=p(0,0)



return
end

!--------------------------------------------------------------
subroutine miscellaneous 
use size
use other
use fdcommon

do i=1,nx-1
if(mod(i,4).eq.0) then
ssx(i)=28d0/45d0*eps
else if(mod(i,2).eq.0) then
ssx(i)=24d0/45d0*eps
else
ssx(i)=64d0/45d0*eps
end if
enddo
ssx(0)=14d0/45d0*eps
ssx(nx)=ssx(0)

do i=1,ny-1
if(mod(i,4).eq.0) then
ssy(i)=28d0/45d0*eps
else if(mod(i,2).eq.0) then
ssy(i)=24d0/45d0*eps
else
ssy(i)=64d0/45d0*eps
end if
enddo
ssy(0)=14d0/45d0*eps
ssy(ny)=ssy(0)

return
end
!--------------------------------------------------------------
subroutine p_shift
use fdcommon
use uvp
use other
implicit none
integer(4)::i,j
real(8)::p_average

p_average=0d0
do j=0,ny
do i=0,nx
p_average=p_average+p(i,j)*ssx(i)*ssy(j)
!write(*,*) i,j,p(i,j)
enddo
enddo
!
p_average=p_average/(2d0*pi)**2

!
do j=0,ny
do i=0,nx
p(i,j)=p(i,j)-p_average
enddo
enddo
return
end


!==========================================
subroutine set_fdcommon
use fdcommon
implicit none
veps=1d0/eps
hveps=veps*0.5d0
vveps=veps**2
hvveps=vveps*.5d0
qvveps=vveps*0.25d0
veps12=veps/12d0
vveps12=vveps/12d0
hvvveps=veps**3*.5d0
vvvveps=veps**4
vbvveps=1d0/beta*vveps
return
end

REAL*8 Function del_x(h,i,j,nx,ny)
use fdcommon
integer,intent(in)::nx,ny
real(8),dimension(-2:nx+1,-2:ny+1),intent(in)::h
!
!   hveps=(1/2)/eps 
!
del_x=(h(i+1,j)-h(i-1,j))*hveps
return
end
!---------------------------------------------------

REAL*8 Function del_xx(h,i,j,nx,ny)
use fdcommon
integer,intent(in)::nx,ny
real(8),dimension(-2:nx+1,-2:ny+1),intent(in)::h
!
!   vveps=1/eps**2 
!
del_xx=(h(i+1,j)-2d0*h(i,j)+h(i-1,j))*vveps
return
end
!----------------------------------------------------

REAL*8 Function del_y(h,i,j,nx,ny)
use fdcommon
integer,intent(in)::nx,ny
real(8),dimension(-2:nx+1,-2:ny+1),intent(in)::h
!
!   hveps=(1/2)/eps 
!
del_y=(h(i,j+1)-h(i,j-1))*hveps
return
end
!------------------------------------------------------

REAL*8 Function del_yy(h,i,j,nx,ny)
use fdcommon
integer,intent(in)::nx,ny
real(8),dimension(-2:nx+1,-2:ny+1),intent(in)::h
!
!   vveps=1/eps**2 
!
del_yy=(h(i,j+1)-2d0*h(i,j)+h(i,j-1))*vveps
return
end
!-------------------------------------------------------

REAL*8 Function del_xy(h,i,j,nx,ny)
use fdcommon
integer,intent(in)::nx,ny
real(8),dimension(-2:nx+1,-2:ny+1),intent(in)::h

!
!   qvveps=(1/4)/eps**2 
!
del_xy=((h(i+1,j+1)-h(i+1,j-1))-(h(i-1,j+1)-h(i-1,j-1)))*qvveps
return
end
!-------------------------------------------------------

REAL*8 Function del_xyy(h,i,j,nx,ny)
use fdcommon
integer,intent(in)::nx,ny
real(8),dimension(-2:nx+1,-2:ny+1),intent(in)::h
!
!   hvvveps=(1/2)/eps**3 
!
del_xyy=((h(i+1,j+1)-2d0*h(i+1,j)+h(i+1,j-1))-(h(i-1,j+1)-2d0*h(i-1,j)+h(i-1,j-1)))*hvvveps
return
end
!----------------------------------------------------------

REAL*8 Function del_xxy(h,i,j,nx,ny)
use fdcommon
integer,intent(in)::nx,ny
real(8),dimension(-2:nx+1,-2:ny+1),intent(in)::h
!
!   hvvveps=(1/2)/eps**3 
!
del_xxy=((h(i+1,j+1)-2d0*h(i,j+1)+h(i-1,j+1))-(h(i+1,j-1)-2d0*h(i,j-1)+h(i-1,j-1)))*hvvveps
return
end
!----------------------------------------------------------

REAL*8 Function del_xxyy(h,i,j,nx,ny)
use fdcommon
integer,intent(in)::nx,ny
real(8),dimension(-2:nx+1,-2:ny+1),intent(in)::h
!
!   vveps=1/eps**4 
!
del_xxyy=((h(i+1,j+1)-2d0*h(i,j+1)+h(i-1,j+1))-2d0*(h(i+1,j)-2d0*h(i,j)+h(i-1,j))+(h(i+1,j-1)-2d0*h(i,j-1)+h(i-1,j-1)))*vvvveps
return
end
!----------------------------------------------------------