|
|
| (22 revisões intermediárias pelo mesmo usuário não estão sendo mostradas) |
| Linha 1: |
Linha 1: |
| <source lang = "fortran">
| | [[Código 4]] |
| MODULE globais
| |
| 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>
| |