Exemplo 4 - Caixa com carga pontual: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
Sem resumo de edição
(Página substituída por '== Código == == Gnuplot ==')
Linha 1: Linha 1:
      MODULE globais
== Código ==
          IMPLICIT NONE
          INTEGER :: N,L
          REAL :: dx
      END MODULE


      PROGRAM exemplo1_poisson
== Gnuplot ==
      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

Edição das 23h19min de 5 de fevereiro de 2023

Código

Gnuplot

Disponível em “http://fiscomp.if.ufrgs.br/index.php?title=Exemplo_4_-_Caixa_com_carga_pontual&oldid=9337”