Código 4: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
Sem resumo de edição
Sem resumo de edição
 
(6 revisões intermediárias pelo mesmo usuário não estão sendo mostradas)
Linha 1: Linha 1:
<source lang = "FORTRAN90">
       MODULE globais
       MODULE globais
           IMPLICIT NONE
           IMPLICIT NONE
Linha 324: Linha 322:
           END DO
           END DO
       END SUBROUTINE dados
       END SUBROUTINE dados
</source>

Edição atual tal como às 23h29min de 5 de fevereiro de 2023

      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