- #1
Nikolas_Ex_Aguirre
- 1
- 0
- TL;DR Summary
- I have a problem with this code that i made. every time that i run it, this shows an error (is the same in the image), and i don't know how to fix this. It's a jacobi relaxation method for a dirichlet boundary conditions.
Code:
program r_jacobi
implicit none
!!!!Variables!!!
real*8 V, V_1, V_2, Lx, Ly
integer n ,i , j, k, nx, ny
real*8, allocatable :: arrx(:), arry(:), phi(:,:,:)
real*8 x, xi, xf, y, yi, yf, dx, dy
real*8 d, q, bx, by
V=1
V_1=V
V_2=-V
Lx = 2
Ly = 1
nx = 200
ny = nx/2
n = 200
k = 200
allocate(arrx(nx))
allocate(arry(ny))
allocate(phi(n,nx,ny))!!!GRID!!!
xi = 0
xf = Lx
dx = (xf-xi)/(nx-1)
yi = 0
yf = Ly
dy = (yf-yi)/(ny-1)
d = 2*(dx**(-1)+dy**(-1)) do j = 1,ny
do i = 1,nx
x = xi+(i-1)*dx
arrx(i) = x
y = yi+(j-1)*dy
arry(j)=y
end do
end do
!!!BC!!!
do j = 1,ny
do i = 1,nx
phi(1,i,j) = (V_1+V_2)/2
end do
end do
do i = 1, nx/2
phi(1,i,1) = V_1
phi(1,i,ny) = V_1
end do
do i = nx/2 +1, nx
phi(1,i,1) = V_2
phi(1,i,ny) = V_2
end do
do j = 1, ny
phi(1,1,j) = V_1
phi(1,nx,j) = V_2
end do
!!!jacobi!!!
!delta_2^x = u^n_j+1_k - 2 u^n_j_k + u^n_j-1_k
!delta_2^y = u^n_j_k+1 - 2 u^n_j_k + u^n_j_k-1
!alpha^0_jk = -2(1/dx^2 + 1/dy^2)
!alpha^1_jk = alpha^2_jk = -1/dx^2
!alpha^3_jk = alpha^4_jk = -1/dy^2 do k = 1, n
do j = 2, ny-1
do i = 2, nx-1
bx = phi(k,i+1,j) - 2*phi(k,i,j) + phi(k,i-1,j)
by = phi(k,i,j+1) - 2*phi(k,i,j) + phi(k,i,j-1)
q = -dx**(-2)*(phi(k,i,j))*bx - dy**(-2)*phi(k,i,j)*by
phi(k+1,i,j)=d**(-1)(q+dx**(-1)*(phi(k,i+1,j)+phi(k,i-1,j))+dy**(-1)*(phi(k,i,j+1)&
+phi(k,i,j-1)))
end do
end do
end do
end program r_jacobi
Last edited by a moderator: