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
Sem resumo de edição
 
(9 revisões intermediárias pelo mesmo usuário não estão sendo mostradas)
Linha 1: Linha 1:
<source lang='fortran'> MODULE globais
[[Código 4]]
          IMPLICIT NONE
          INTEGER :: N,L
          REAL :: dx
      END MODULE


      PROGRAM exemplo1_poisson
[[Gnuplot 4]]
      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
</source>

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

Código 4

Gnuplot 4

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