Código 4
Ir para navegação
Ir para pesquisar
MODULE globais
IMPLICIT NONE
INTEGER :: N,L
REAL :: dx
END MODULE
PROGRAM exemplo1_poisson
USE globais
IMPLICIT NONE
INTEGER :: i,j,k,c,t,edge
REAL :: dV
REAL,DIMENSION(:,:,:),ALLOCATABLE :: V,V1,E1,E2,E3,x,y,z,dens
N = 1
edge = 20*N
L = edge + 1
dx = 0.1
ALLOCATE(V(L,L,L),V1(L,L,L),E1(L,L,L),E2(L,L,L),E3(L,L,L), &
x(L,L,L),y(L,L,L),z(L,L,L),dens(L,L,L))
DO k=1,L
DO i=1,L
DO j=1,L
x(j,i,k) = -1.0 + (i-1)*dx/N
y(j,i,k) = -1.0 + (j-1)*dx/N
z(j,i,k) = -1.0 + (k-1)*dx/N
END DO
END DO
END DO
! Calcula utilizando o algoritmo de Jacobi
CALL inicia_zeros(E1)
CALL inicia_zeros(E2)
CALL inicia_zeros(E3)
CALL initialize_V(V,dens)
c = 0
V1 = V
CALL laplace_jacobi(V,V1,dens,dV,c)
CALL electric_field(V,E1,E2,E3)
WRITE(*,*) 'Jacobi:',dV,c
CALL dados(x,y,z,V,E1,E2,E3)
! Calcula utilizando o algoritmo de Gauss-Siedel
CALL inicia_zeros(E1)
CALL inicia_zeros(E2)
CALL inicia_zeros(E3)
CALL initialize_V(V,dens)
c = 0
CALL laplace_gauss_siedel(V,dens,dV,c)
CALL electric_field(V,E1,E2,E3)
WRITE(*,*) 'Gauss-Siedel:',dV,c
CALL dados(x,y,z,V,E1,E2,E3)
END PROGRAM
SUBROUTINE inicia_zeros(M)
! Inicializa uma matriz em zeros
USE globais
IMPLICIT NONE
INTEGER :: i,j,k
REAL,DIMENSION(L,L,L) :: M
DO k=1,L
DO i=1,L
DO j=1,L
M(j,i,k) = 0.0
END DO
END DO
END DO
END SUBROUTINE inicia_zeros
SUBROUTINE initialize_V(V0,rho)
! Inicialização do potencial
USE globais
IMPLICIT NONE
INTEGER :: i,j,k
REAL,DIMENSION(L,L,L) :: V0,rho
! Inicializa o potencial
DO k=1,L
DO i=1,L
DO j=1,L
V0(j,i,k) = 0.0
rho(j,i,k) = 0.0
END DO
END DO
END DO
rho(11,11,11) = 1.0/dx**3
rho(11,11,10) = 1.0/dx**3
rho(11,11,12) = 1.0/dx**3
rho(11,10,11) = 1.0/dx**3
rho(11,12,11) = 1.0/dx**3
rho(10,11,11) = 1.0/dx**3
rho(12,11,11) = 1.0/dx**3
END SUBROUTINE initialize_V
SUBROUTINE update_V_jacobi(Vn,Vn1,rho,deltaV)
! Calcula a atualização uma vez
USE globais
IMPLICIT NONE
INTEGER :: i,j,k
REAL :: deltaV
REAL,DIMENSION(L,L,L) :: Vn,Vn1,rho
deltaV = 0.0
DO k=2,L-1
DO i=2,L-1
DO j=2,L-1
Vn1(j,i,k) = (Vn(j+1,i,k) + Vn(j-1,i,k) + &
Vn(j,i+1,k) + Vn(j,i-1,k) + Vn(j,i,k+1) + Vn(j,i,k-1))/6.0 &
+ (rho(j,i,k)*dx**2)/6.0
deltaV = deltaV + abs(Vn(j,i,k) - Vn1(j,i,k))
END DO
END DO
END DO
END SUBROUTINE update_V_jacobi
SUBROUTINE laplace_jacobi(Vm,Vm1,rho,delta_v,cont)
! Implementa o método de Jacobi
USE globais
IMPLICIT NONE
INTEGER :: cont
REAL :: delta_v,tol
REAL,DIMENSION(L,L,L) :: Vm,Vm1,rho
tol = 1e-5*L*L*L
delta_v = 100.0
DO WHILE (delta_v .gt. tol)
CALL update_V_jacobi(Vm,Vm1,rho,delta_v)
CALL update_V_jacobi(Vm1,Vm,rho,delta_v)
cont = cont+1
END DO
END SUBROUTINE laplace_jacobi
SUBROUTINE update_V_gauss_siedel(Vn,rho,deltaV)
! Calcula a atualização uma vez
USE globais
IMPLICIT NONE
INTEGER :: i,j,k
REAL :: deltaV,v1,v2
REAL,DIMENSION(L,L,L) :: Vn,rho
deltaV = 0.0
DO k=2,L-1
DO i=2,L-1
DO j=2,L-1
v1 = Vn(j,i,k)
Vn(j,i,k) = (Vn(j+1,i,k) + Vn(j-1,i,k) + &
Vn(j,i+1,k) + Vn(j,i-1,k) + Vn(j,i,k+1) + Vn(j,i,k-1))/6.0 &
+ (rho(j,i,k)*dx**2)/6.0
v2 = Vn(j,i,k)
deltaV = deltaV + abs(v2 - v1)
END DO
END DO
END DO
END SUBROUTINE update_V_gauss_siedel
SUBROUTINE laplace_gauss_siedel(Vm,rho,delta_v,cont)
! Implementa o método de Gauss-Siedel
USE globais
IMPLICIT NONE
INTEGER :: cont
REAL :: delta_v,tol
REAL,DIMENSION(L,L,L) :: Vm,rho
tol = 1e-5*L*L*L
delta_v = 100.0
DO WHILE (delta_v .gt. tol)
CALL update_V_gauss_siedel(Vm,rho,delta_v)
CALL update_V_gauss_siedel(Vm,rho,delta_v)
cont = cont+1
END DO
END SUBROUTINE laplace_gauss_siedel
SUBROUTINE electric_field(Vk,Ex,Ey,Ez)
! Calcula o campo elétrico a partir do potencial
USE globais
IMPLICIT NONE
INTEGER :: i,j,k
REAL,DIMENSION(L,L,L) :: Vk,Ex,Ey,Ez
! Bordas k=1 e k=L
DO k=1,L,L-1
DO i=1,L
DO j=1,L
! Derivada para trás
IF (i == L) THEN
Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k))
Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k))
Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1))
! Derivada para trás
ELSE IF (j == L) THEN
Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k))
Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k))
Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1))
! Derivada para trás
ELSE IF (k == L) THEN
Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k))
Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k))
Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1))
! Derivada para frente
ELSE
Ex(j,i,k) = - (Vk(j,i+1,k) - Vk(j,i,k))
Ey(j,i,k) = - (Vk(j+1,i,k) - Vk(j,i,k))
Ez(j,i,k) = - (Vk(j,i,k+1) - Vk(j,i,k))
END IF
END DO
END DO
END DO
! Bordas j=1 e j=L
DO j=1,L,L-1
DO i=1,L
DO k=1,L
! Derivada para trás
IF (i == L) THEN
Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k))
Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k))
Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1))
! Derivada para trás
ELSE IF (j == L) THEN
Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k))
Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k))
Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1))
! Derivada para trás
ELSE IF (k == L) THEN
Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k))
Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k))
Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1))
! Derivada para frente
ELSE
Ex(j,i,k) = - (Vk(j,i+1,k) - Vk(j,i,k))
Ey(j,i,k) = - (Vk(j+1,i,k) - Vk(j,i,k))
Ez(j,i,k) = - (Vk(j,i,k+1) - Vk(j,i,k))
END IF
END DO
END DO
END DO
! Bordas i=1 e i=L
DO i=1,L,L-1
DO k=1,L
DO j=1,L
! Derivada para trás
IF (i == L) THEN
Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k))
Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k))
Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1))
! Derivada para trás
ELSE IF (j == L) THEN
Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k))
Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k))
Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1))
! Derivada para trás
ELSE IF (k == L) THEN
Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k))
Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k))
Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1))
! Derivada para frente
ELSE
Ex(j,i,k) = - (Vk(j,i+1,k) - Vk(j,i,k))
Ey(j,i,k) = - (Vk(j+1,i,k) - Vk(j,i,k))
Ez(j,i,k) = - (Vk(j,i,k+1) - Vk(j,i,k))
END IF
END DO
END DO
END DO
! Derivada centrada
DO k=2,L-1
DO i=2,L-1
DO j=2,L-1
Ex(j,i,k) = - (Vk(j,i+1,k) - Vk(j,i-1,k))/2.0
Ey(j,i,k) = - (Vk(j+1,i,k) - Vk(j-1,i,k))/2.0
Ez(j,i,k) = - (Vk(j,i,k+1) - Vk(j,i,k-1))/2.0
END DO
END DO
END DO
END SUBROUTINE electric_field
SUBROUTINE dados(x1,y1,z1,Vk,Ex,Ey,Ez)
! Imprime os dados em arquivos
USE globais
IMPLICIT NONE
INTEGER :: i,j,k,m
REAL,DIMENSION(L,L,L) :: Vk,Ex,Ey,Ez,x1,y1,z1
CHARACTER(28),DIMENSION(2) :: a
CHARACTER(27),DIMENSION(2) :: b
a(1) = 'potencial_poisson_jacob1.dat'
a(2) = 'potencial_poisson_gauss1.dat'
b(1) = 'electric_poisson_jacob1.dat'
b(2) = 'electric_poisson_gauss1.dat'
DO m=1,2
OPEN(10,FILE=a(m),STATUS='unknown',ACTION='write')
DO i=1,L
DO j=1,L
DO k=1,L
WRITE(10,*) x1(j,i,k),y1(j,i,k),z1(j,i,k),Vk(j,i,k)
END DO
WRITE(10,*)
END DO
WRITE(10,*)
END DO
CLOSE(10)
OPEN(20,FILE=b(m),STATUS='unknown',ACTION='write')
DO i=1,L
DO j=1,L
DO k=1,L
WRITE(20,*) x1(j,i,k),y1(j,i,k),z1(j,i,k), &
Ex(j,i,k),Ey(j,i,k),Ez(j,i,k)
END DO
WRITE(20,*)
END DO
WRITE(20,*)
END DO
CLOSE(20)
END DO
END SUBROUTINE dados