Exibir código-fonte para 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