Exemplo 4 - Caixa com carga pontual
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