program wave_single
implicit none
integer , parameter :: nx = 1024
integer , parameter :: ny = 1024
integer , parameter :: b_nx = 8
integer , parameter :: b_ny = 8
integer , parameter :: itterations = 10
integer :: i , j, k, itt, max_itt
double precision :: dx , dy, dt
double precision :: xmin , xmax, ymin, ymax
double precision :: r , t_start, t_end
double precision :: res , jac
real :: start_time , stop_time
double precision , dimension ( nx, ny) :: y0_n , y1_n, y0_np1
double precision , dimension ( nx,ny) :: y1_np1 , temp0, temp1
double precision , dimension ( nx) :: x
double precision , dimension ( ny) :: y
xmin = - 10.0
xmax = 10.0
ymin = - 10.0
ymax = 10.0
max_itt = itterations * nx/ b_nx
dx = ( xmax- xmin) / DBLE ( nx- 1 )
dy = ( ymax- ymin) / DBLE ( ny- 1 )
dt = 0.1 * ( xmax- xmin) / DBLE ( nx)
do i = 1 , nx
x( i) = xmin + ( i- 1 ) * dx
enddo
do j = 1 , ny
y( j) = ymin + ( j- 1 ) * dy
enddo
do i = 1 , nx
do j = 1 , ny
y0_n( i,j) = 0
y1_n( i,j) = 0
y0_np1( i,j) = 0
y1_np1( i,j) = 0
enddo
enddo
do i = 2 , nx- 1
do j = 2 , ny- 1
y0_n( i,j) = dexp( - x( i) * x( i) - y( j) * y( j) )
y0_np1( i,j) = dexp( - x( i) * x( i) - y( j) * y( j) )
enddo
enddo
call cpu_time ( start_time)
do itt = 0 , max_itt - 1
do k = 0 , 10
do i = 2 , nx- 1
do j = 2 , ny- 1
res = ( y0_np1( i,j) - y0_n( i,j) ) / dt &
- 0.5 * ( y1_np1( i,j) + y1_n( i,j) )
jac = 1d0/ dt
temp0( i,j) = y0_np1( i,j) - res/ jac
res = ( y1_np1( i,j) - y1_n( i,j) ) / dt &
- 0.5 * ( y0_np1( i- 1 ,j) - 2 * y0_np1( i,j) + y0_np1( i+ 1 ,j) ) / ( dx** 2 ) &
- 0.5 * ( y0_np1( i,j- 1 ) - 2 * y0_np1( i,j) + y0_np1( i,j+ 1 ) ) / ( dy** 2 ) &
- 0.5 * ( y0_n( i- 1 ,j) - 2 * y0_n( i,j) + y0_n( i+ 1 ,j) ) / ( dx** 2 ) &
- 0.5 * ( y0_n( i,j- 1 ) - 2 * y0_n( i,j) + y0_n( i,j+ 1 ) ) / ( dy** 2 )
jac = 1d0/ dt
temp1( i,j) = y1_np1( i,j) - res/ jac
enddo
enddo
do i = 2 , nx- 1
do j = 2 , ny- 1
y0_np1( i,j) = temp0( i,j)
y1_np1( i,j) = temp1( i,j)
enddo
enddo
enddo
do i = 1 , nx
do j = 1 , ny
res = y0_n( i,j)
y0_n( i,j) = y0_np1( i,j)
y0_np1( i,j) = res
res = y1_n( i,j)
y1_n( i,j) = y1_np1( i,j)
y1_np1( i,j) = res
enddo
enddo
enddo
call cpu_time ( stop_time)
write( * ,* ) "total time: " , stop_time- start_time
!call print2d(y0_n,nx,ny)
!call print2d(y1_n,nx,ny)
!call print2d(y0_np1,nx,ny)
!call print2d(y1_np1,nx,ny)
end
subroutine print2d( y, nx, ny)
implicit none
integer :: nx , ny, i, j
double precision , dimension ( nx,ny) :: y
do j = 1 , ny
do i = 1 , nx
if ( y( i,j) > 0.05 ) then
write( * ,"(i2)" ,advance = "no" ) 1
elseif ( y( i,j) < - 0.05 ) then
write( * ,"(i2)" ,advance = "no" ) - 1
else
write( * ,"(i2)" , advance = "no" ) 0
endif
enddo
write( * ,* ) ""
enddo
write( * ,* ) ""
end
cHJvZ3JhbSB3YXZlX3NpbmdsZQogICAgICBpbXBsaWNpdCBub25lCiAgICAgIGludGVnZXIsIHBhcmFtZXRlciAgICAgICA6OiBueCA9IDEwMjQKICAgICAgaW50ZWdlciwgcGFyYW1ldGVyICAgICAgIDo6IG55ID0gMTAyNAogICAgICBpbnRlZ2VyLCBwYXJhbWV0ZXIgICAgICAgOjogYl9ueCA9IDgKICAgICAgaW50ZWdlciwgcGFyYW1ldGVyICAgICAgIDo6IGJfbnkgPSA4CiAgICAgIGludGVnZXIsIHBhcmFtZXRlciAgICAgICA6OiBpdHRlcmF0aW9ucyA9IDEwCgogICAgICBpbnRlZ2VyICAgICAgICAgICAgICAgICAgOjogaSwgaiwgaywgaXR0LCBtYXhfaXR0CiAgICAgIGRvdWJsZSBwcmVjaXNpb24gICAgICAgICA6OiBkeCwgZHksIGR0CiAgICAgIGRvdWJsZSBwcmVjaXNpb24gICAgICAgICA6OiB4bWluLCB4bWF4LCB5bWluLCB5bWF4CiAgICAgIGRvdWJsZSBwcmVjaXNpb24gICAgICAgICA6OiByLCB0X3N0YXJ0LCB0X2VuZAogICAgICBkb3VibGUgcHJlY2lzaW9uICAgICAgICAgOjogcmVzLCBqYWMKICAgICAgcmVhbCAgICAgICAgICAgICAgICAgICAgIDo6IHN0YXJ0X3RpbWUsIHN0b3BfdGltZQoKICAgICAgZG91YmxlIHByZWNpc2lvbiwgZGltZW5zaW9uKG54LCBueSkgOjogeTBfbiwgeTFfbiwgeTBfbnAxCiAgICAgIGRvdWJsZSBwcmVjaXNpb24sIGRpbWVuc2lvbihueCxueSkgIDo6IHkxX25wMSwgdGVtcDAsIHRlbXAxCiAgICAgIGRvdWJsZSBwcmVjaXNpb24sIGRpbWVuc2lvbihueCkgICAgIDo6IHgKICAgICAgZG91YmxlIHByZWNpc2lvbiwgZGltZW5zaW9uKG55KSAgICAgOjogeQoKCiAgICAgIHhtaW4gPSAtMTAuMAogICAgICB4bWF4ID0gMTAuMAogICAgICB5bWluID0gLTEwLjAKICAgICAgeW1heCA9IDEwLjAKICAgICAgbWF4X2l0dCA9IGl0dGVyYXRpb25zICogbngvYl9ueAoKICAgICAgZHggPSAoeG1heC14bWluKS9EQkxFKG54LTEpCiAgICAgIGR5ID0gKHltYXgteW1pbikvREJMRShueS0xKQogICAgICBkdCA9IDAuMSooeG1heC14bWluKS9EQkxFKG54KQogICAgICAKICAgICAgZG8gaSA9IDEsIG54CiAgICAgICB4KGkpID0geG1pbiArIChpLTEpKmR4CiAgICAgIGVuZGRvCgogICAgICBkbyBqID0gMSwgbnkKICAgICAgIHkoaikgPSB5bWluICsgKGotMSkqZHkKICAgICAgZW5kZG8KCiAgICAgIGRvIGkgPSAxLCBueAogICAgICAgZG8gaiA9IDEsIG55CiAgICAgICAgeTBfbihpLGopICAgPSAwCiAgICAgICAgeTFfbihpLGopICAgPSAwCiAgICAgICAgeTBfbnAxKGksaikgPSAwCiAgICAgICAgeTFfbnAxKGksaikgPSAwIAogICAgICAgZW5kZG8KICAgICAgZW5kZG8KCiAgICAgIGRvIGkgPSAyLCBueC0xCiAgICAgICBkbyBqID0gMiwgbnktMQogICAgICAgIHkwX24oaSxqKSAgID0gZGV4cCgteChpKSp4KGkpLXkoaikqeShqKSkKICAgICAgICB5MF9ucDEoaSxqKSA9IGRleHAoLXgoaSkqeChpKS15KGopKnkoaikpCiAgICAgICBlbmRkbwogICAgICBlbmRkbyAgICAgIAoKICAgICAgY2FsbCBjcHVfdGltZShzdGFydF90aW1lKSAKCiAgICAgIGRvIGl0dCA9IDAsIG1heF9pdHQgLSAxCiAgICAgICBkbyBrID0gMCwgMTAKICAgICAgICBkbyBpID0gMiwgbngtMQogICAgICAgICBkbyBqID0gMiwgbnktMQogICAgICAgICAgcmVzID0gKHkwX25wMShpLGopIC0geTBfbihpLGopKS9kdCAmIAogICAgICAgLSAwLjUqKHkxX25wMShpLGopK3kxX24oaSxqKSkKICAgICAgICAgIGphYyA9IDFkMC9kdAogICAgICAgICAgdGVtcDAoaSxqKSA9IHkwX25wMShpLGopIC0gcmVzL2phYwoKICAgICAgICAgIHJlcyA9ICh5MV9ucDEoaSxqKSAtIHkxX24oaSxqKSkvZHQgJiAKICAgICAgIC0wLjUqKHkwX25wMShpLTEsaiktMip5MF9ucDEoaSxqKSt5MF9ucDEoaSsxLGopKS8oZHgqKjIpICYgCiAgICAgICAtMC41Kih5MF9ucDEoaSxqLTEpLTIqeTBfbnAxKGksaikreTBfbnAxKGksaisxKSkvKGR5KioyKSAmIAogICAgICAgLTAuNSooeTBfbihpLTEsaiktMip5MF9uKGksaikreTBfbihpKzEsaikpLyhkeCoqMikgJiAKICAgICAgIC0wLjUqKHkwX24oaSxqLTEpLTIqeTBfbihpLGopK3kwX24oaSxqKzEpKS8oZHkqKjIpCiAgICAgICAgICBqYWMgPSAxZDAvZHQKICAgICAgICAgIHRlbXAxKGksaikgPSB5MV9ucDEoaSxqKSAtIHJlcy9qYWMKICAgICAgICAgZW5kZG8KICAgICAgICBlbmRkbwoKICAgICAgICBkbyBpID0gMiwgbngtMQogICAgICAgICBkbyBqID0gMiwgbnktMQogICAgICAgICAgeTBfbnAxKGksaikgPSB0ZW1wMChpLGopCiAgICAgICAgICB5MV9ucDEoaSxqKSA9IHRlbXAxKGksaikKICAgICAgICAgZW5kZG8KICAgICAgICBlbmRkbwogICAgICAgZW5kZG8KICAgICAgIAogICAgICAgZG8gaSA9IDEsIG54CiAgICAgICAgZG8gaiA9IDEsIG55CiAgICAgICAgIHJlcyAgICAgICAgID0geTBfbihpLGopCiAgICAgICAgIHkwX24oaSxqKSAgID0geTBfbnAxKGksaikKICAgICAgICAgeTBfbnAxKGksaikgPSByZXMKCiAgICAgICAgIHJlcyAgICAgICAgID0geTFfbihpLGopCiAgICAgICAgIHkxX24oaSxqKSAgID0geTFfbnAxKGksaikKICAgICAgICAgeTFfbnAxKGksaikgPSByZXMKICAgICAgICBlbmRkbwogICAgICAgZW5kZG8KICAgICAgZW5kZG8KCiAgICAgIGNhbGwgY3B1X3RpbWUoc3RvcF90aW1lKSAKICAgICAgd3JpdGUoKiwqKSAidG90YWwgdGltZTogIiwgc3RvcF90aW1lLXN0YXJ0X3RpbWUKCiAgICAgICFjYWxsIHByaW50MmQoeTBfbixueCxueSkKICAgICAgIWNhbGwgcHJpbnQyZCh5MV9uLG54LG55KQogICAgICAhY2FsbCBwcmludDJkKHkwX25wMSxueCxueSkKICAgICAgIWNhbGwgcHJpbnQyZCh5MV9ucDEsbngsbnkpCiAgICAgIAogICAgICBlbmQKCiAgICAgIHN1YnJvdXRpbmUgcHJpbnQyZCh5LCBueCwgbnkpCiAgICAgIGltcGxpY2l0IG5vbmUKICAgICAgaW50ZWdlciAgICAgICAgICAgICAgICAgICAgICAgICAgICA6OiAgbngsIG55LCBpLCBqCiAgICAgIGRvdWJsZSBwcmVjaXNpb24sIGRpbWVuc2lvbihueCxueSkgOjogIHkKCiAgICAgIGRvIGogPSAxLCBueQogICAgICAgZG8gaSA9IDEsIG54CiAgICAgICAgaWYgKHkoaSxqKSA+IDAuMDUpIHRoZW4gICAKICAgICAgICAgd3JpdGUoKiwiKGkyKSIsYWR2YW5jZT0ibm8iKSAxCiAgICAgICAgZWxzZWlmICh5KGksaikgPCAtMC4wNSkgdGhlbgogICAgICAgICB3cml0ZSgqLCIoaTIpIixhZHZhbmNlPSJubyIpIC0xCiAgICAgICAgZWxzZQogICAgICAgICB3cml0ZSgqLCIoaTIpIiwgYWR2YW5jZT0ibm8iKSAwCiAgICAgICAgZW5kaWYKICAgICAgIGVuZGRvCiAgICAgICB3cml0ZSgqLCopICIiCiAgICAgIGVuZGRvCiAgICAgIHdyaXRlKCosKikgIiIKCiAgICAgIGVuZA==