Código 4

De Física Computacional
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