Código 3
Ir para navegação
Ir para pesquisar
MODULE globais
IMPLICIT NONE
INTEGER :: N,L
END MODULE
PROGRAM exemplo3_laplace
USE globais
IMPLICIT NONE
INTEGER :: i,j,k,c,t,edge
REAL :: dV
INTEGER,DIMENSION(:,:),ALLOCATABLE :: masc
REAL,DIMENSION(:,:),ALLOCATABLE :: V,V1,E1,E2,x,y
N = 10
edge = 20*N
L = edge + 1
ALLOCATE(V(L,L),V1(L,L),E1(L,L),E2(L,L),x(L,L),y(L,L),masc(L,L))
DO i=1,L
DO j=1,L
x(j,i) = -1.0 + (i-1)*0.1/N
y(j,i) = -1.0 + (j-1)*0.1/N
END DO
END DO
! Calcula utilizando o algoritmo de Jacobi
CALL inicia_zeros(E1)
CALL inicia_zeros(E2)
CALL initialize_V(V,masc)
c = 0
V1 = V
CALL laplace_jacobi(V,V1,masc,dV,c)
CALL electric_field(V,E1,E2)
WRITE(*,*) 'Jacobi:',dV,c
CALL dados(x,y,V,E1,E2)
! Calcula utilizando o algoritmo de Gauss-Siedel
CALL inicia_zeros(E1)
CALL inicia_zeros(E2)
CALL initialize_V(V,masc)
c = 0
CALL laplace_gauss_siedel(V,masc,dV,c)
CALL electric_field(V,E1,E2)
WRITE(*,*) 'Gauss-Siedel:',dV,c
CALL dados(x,y,V,E1,E2)
END PROGRAM
SUBROUTINE inicia_zeros(M)
! Inicializa uma matriz em zeros
USE globais
IMPLICIT NONE
INTEGER :: i,j
REAL,DIMENSION(L,L) :: M
DO i=1,L
DO j=1,L
M(j,i) = 0.0
END DO
END DO
END SUBROUTINE inicia_zeros
SUBROUTINE initialize_V(V0,mascara)
! Inicialização do potencial
USE globais
IMPLICIT NONE
INTEGER :: i,j,m,o,p,q
INTEGER,DIMENSION(L,L) :: mascara
REAL,DIMENSION(L,L) :: V0
! Inicializa a mascara em 0
DO m=1,L
DO o=1,L
mascara(o,m) = 0
END DO
END DO
! Escreve 1's e -1's onde o V=1.0 ou V=-1.0
DO q=7*N+1,13*N+1
mascara(q,71) = 1
mascara(q,131) = -1
END DO
! Escreve o potencial
DO i=1,L
DO j=1,L
IF (mascara(j,i) == 1) THEN
V0(j,i) = 1.0
ELSE IF (mascara(j,i) == -1) THEN
V0(j,i) = -1.0
ELSE
V0(j,i) = 0.0
END IF
END DO
END DO
END SUBROUTINE initialize_V
SUBROUTINE update_V_jacobi(Vn,Vn1,mascara,deltaV)
! Calcula a atualização uma vez
USE globais
IMPLICIT NONE
INTEGER :: i,j
REAL :: deltaV
INTEGER,DIMENSION(L,L) :: mascara
REAL,DIMENSION(L,L) :: Vn,Vn1
deltaV = 0.0
DO i=2,L-1
DO j=2,L-1
IF (mascara(j,i) == 0) THEN
Vn1(j,i) = (Vn(j+1,i) + Vn(j-1,i) + Vn(j,i+1) + &
Vn(j,i-1))/4.0
deltaV = deltaV + abs(Vn(j,i) - Vn1(j,i))
END IF
END DO
END DO
END SUBROUTINE update_V_jacobi
SUBROUTINE laplace_jacobi(Vm,Vm1,mascara,delta_v,cont)
! Implementa o método de Jacobi
USE globais
IMPLICIT NONE
INTEGER :: cont,x
REAL :: delta_v,tol
INTEGER,DIMENSION(L,L) :: mascara
REAL,DIMENSION(L,L) :: Vm,Vm1
tol = 1e-5*L*L
delta_v = 41.0
DO WHILE (delta_v .gt. tol)
CALL update_V_jacobi(Vm,Vm1,mascara,delta_v)
CALL update_V_jacobi(Vm1,Vm,mascara,delta_v)
cont = cont+1
END DO
END SUBROUTINE laplace_jacobi
SUBROUTINE update_V_gauss_siedel(Vn,mascara,deltaV)
! Calcula a atualização uma vez
USE globais
IMPLICIT NONE
INTEGER :: i,j
REAL :: deltaV,v1,v2
INTEGER,DIMENSION(L,L) :: mascara
REAL,DIMENSION(L,L) :: Vn
deltaV = 0.0
DO i=2,L-1
DO j=2,L-1
IF (mascara(j,i) == 0) THEN
v1 = Vn(j,i)
Vn(j,i) = (Vn(j+1,i) + Vn(j-1,i) + Vn(j,i+1) + &
Vn(j,i-1))/4.0
v2 = Vn(j,i)
deltaV = deltaV + abs(v2 - v1)
END IF
END DO
END DO
END SUBROUTINE update_V_gauss_siedel
SUBROUTINE laplace_gauss_siedel(Vm,mascara,delta_v,cont)
! Implementa o método de Gauss_Siedel
USE globais
IMPLICIT NONE
INTEGER :: cont
REAL :: delta_v,tol
INTEGER,DIMENSION(L,L) :: mascara
REAL,DIMENSION(L,L) :: Vm
tol = 1e-5*L*L
delta_v = 10.0
DO WHILE (delta_v .gt. tol)
CALL update_V_gauss_siedel(Vm,mascara,delta_v)
CALL update_V_gauss_siedel(Vm,mascara,delta_v)
cont = cont+1
END DO
END SUBROUTINE laplace_gauss_siedel
SUBROUTINE electric_field(Vk,Ex,Ey)
! Calcula o campo elétrico a partir do potencial
USE globais
IMPLICIT NONE
INTEGER :: i,j
REAL,DIMENSION(L,L) :: Vk,Ex,Ey
! Derivada para frente (bordas)
DO i=1,L-1
Ex(1,i) = - (Vk(1,i+1) - Vk(1,i))
Ex(L,i) = - (Vk(L,i+1) - Vk(L,i))
Ey(i,1) = - (Vk(i+1,1) - Vk(i,1))
Ey(i,L) = - (Vk(i+1,L) - Vk(i,L))
END DO
! Derivada para trás (bordas)
Ex(1,L) = - (Vk(1,L) - Vk(1,L-1))
Ex(L,L) = - (Vk(L,L) - Vk(L,L-1))
Ey(L,1) = - (Vk(L,1) - Vk(L-1,1))
Ey(L,L) = - (Vk(L,L) - Vk(L-1,L))
! Derivada centrada
DO i=2,L-1
DO j=2,L-1
Ex(j,i) = - (Vk(j,i+1) - Vk(j,i-1))/2.0
Ey(j,i) = - (Vk(j+1,i) - Vk(j-1,i))/2.0
END DO
END DO
END SUBROUTINE electric_field
SUBROUTINE dados(x1,y1,Vk,Ex,Ey)
! Imprime os dados em arquivos
USE globais
IMPLICIT NONE
INTEGER :: i,j,k
REAL,DIMENSION(L,L) :: Vk,Ex,Ey,x1,y1
CHARACTER(28),DIMENSION(2) :: a
CHARACTER(27),DIMENSION(2) :: b
CHARACTER(32),DIMENSION(2) :: c
a(1) = 'potencial_laplace_jacob3.dat'
a(2) = 'potencial_laplace_gauss3.dat'
b(1) = 'electric_laplace_jacob3.dat'
b(2) = 'electric_laplace_gauss3.dat'
c(1) = 'equipotencial_laplace_jacob3.dat'
c(2) = 'equipotencial_laplace_gauss3.dat'
DO k=1,2
OPEN(10,FILE=a(k),STATUS='unknown',ACTION='write')
DO i=1,L
DO j=1,L
WRITE(10,*) x1(j,i),y1(j,i),Vk(j,i)
END DO
WRITE(10,*)
END DO
CLOSE(10)
OPEN(20,FILE=b(k),STATUS='unknown',ACTION='write')
DO i=1,L
DO j=1,L
WRITE(20,*) x1(j,i),y1(j,i),Ex(j,i),Ey(j,i)
END DO
WRITE(20,*)
END DO
CLOSE(20)
OPEN(30,FILE=c(k),STATUS='unknown',ACTION='write')
DO i=1,L
DO j=1,L
IF ((Vk(j,i) .lt. 0.96) .and. (Vk(j,i) .gt. 0.94 &
)) THEN
WRITE(30,*) x1(j,i),y1(j,i)
ELSE IF ((Vk(j,i) .lt. 0.76) .and. (Vk(j,i) .gt. &
0.74 )) THEN
WRITE(30,*) x1(j,i),y1(j,i)
ELSE IF ((Vk(j,i) .lt. 0.56) .and. (Vk(j,i) .gt. &
0.54 )) THEN
WRITE(30,*) x1(j,i),y1(j,i)
ELSE IF ((Vk(j,i) .lt. 0.36) .and. (Vk(j,i) .gt. &
0.34 )) THEN
WRITE(30,*) x1(j,i),y1(j,i)
ELSE IF ((Vk(j,i) .lt. 0.16) .and. (Vk(j,i) .gt. &
0.14 )) THEN
WRITE(30,*) x1(j,i),y1(j,i)
ELSE IF ((Vk(j,i) .gt. -0.96) .and. (Vk(j,i) .lt.&
-0.94 )) THEN
WRITE(30,*) x1(j,i),y1(j,i)
ELSE IF ((Vk(j,i) .gt. -0.76) .and. (Vk(j,i) .lt.&
-0.74 )) THEN
WRITE(30,*) x1(j,i),y1(j,i)
ELSE IF ((Vk(j,i) .gt. -0.56) .and. (Vk(j,i) .lt.&
-0.54 )) THEN
WRITE(30,*) x1(j,i),y1(j,i)
ELSE IF ((Vk(j,i) .gt. -0.36) .and. (Vk(j,i) .lt.&
-0.34 )) THEN
WRITE(30,*) x1(j,i),y1(j,i)
ELSE IF ((Vk(j,i) .gt. -0.16) .and. (Vk(j,i) .lt.&
-0.14 )) THEN
WRITE(30,*) x1(j,i),y1(j,i)
END IF
END DO
END DO
CLOSE(30)
END DO
END SUBROUTINE dados