http://fiscomp.if.ufrgs.br/index.php?action=history&feed=atom&title=C%C3%B3digo 2026-04-16T02:49:43Z Histórico de revisões para esta página neste wiki MediaWiki 1.39.4 http://fiscomp.if.ufrgs.br/index.php?title=C%C3%B3digo&diff=9339&oldid=prev 2023-02-06T02:19:59Z

<p>Criou página com &#039; MODULE globais IMPLICIT NONE INTEGER :: N,L REAL :: dx END MODULE PROGRAM exemplo1_poisson USE globais IM...&#039;</p> <p><b>Página nova</b></p><div> MODULE globais<br /> IMPLICIT NONE<br /> INTEGER :: N,L<br /> REAL :: dx <br /> END MODULE<br /> <br /> PROGRAM exemplo1_poisson<br /> USE globais<br /> IMPLICIT NONE<br /> INTEGER :: i,j,k,c,t,edge<br /> REAL :: dV<br /> REAL,DIMENSION(:,:,:),ALLOCATABLE :: V,V1,E1,E2,E3,x,y,z,dens <br /> <br /> N = 1<br /> edge = 20*N<br /> L = edge + 1<br /> dx = 0.1<br /> <br /> ALLOCATE(V(L,L,L),V1(L,L,L),E1(L,L,L),E2(L,L,L),E3(L,L,L), &amp;<br /> x(L,L,L),y(L,L,L),z(L,L,L),dens(L,L,L)) <br /> <br /> DO k=1,L<br /> DO i=1,L<br /> DO j=1,L<br /> x(j,i,k) = -1.0 + (i-1)*dx/N<br /> y(j,i,k) = -1.0 + (j-1)*dx/N<br /> z(j,i,k) = -1.0 + (k-1)*dx/N<br /> END DO<br /> END DO<br /> END DO<br /> <br /> ! Calcula utilizando o algoritmo de Jacobi <br /> CALL inicia_zeros(E1)<br /> CALL inicia_zeros(E2)<br /> CALL inicia_zeros(E3)<br /> CALL initialize_V(V,dens)<br /> <br /> c = 0<br /> V1 = V<br /> <br /> CALL laplace_jacobi(V,V1,dens,dV,c)<br /> CALL electric_field(V,E1,E2,E3)<br /> <br /> WRITE(*,*) 'Jacobi:',dV,c<br /> <br /> CALL dados(x,y,z,V,E1,E2,E3)<br /> <br /> ! Calcula utilizando o algoritmo de Gauss-Siedel<br /> CALL inicia_zeros(E1)<br /> CALL inicia_zeros(E2)<br /> CALL inicia_zeros(E3)<br /> CALL initialize_V(V,dens)<br /> <br /> c = 0<br /> <br /> CALL laplace_gauss_siedel(V,dens,dV,c)<br /> CALL electric_field(V,E1,E2,E3)<br /> <br /> WRITE(*,*) 'Gauss-Siedel:',dV,c<br /> <br /> CALL dados(x,y,z,V,E1,E2,E3)<br /> <br /> END PROGRAM<br /> <br /> SUBROUTINE inicia_zeros(M)<br /> ! Inicializa uma matriz em zeros<br /> USE globais<br /> IMPLICIT NONE<br /> INTEGER :: i,j,k<br /> REAL,DIMENSION(L,L,L) :: M<br /> DO k=1,L<br /> DO i=1,L<br /> DO j=1,L<br /> M(j,i,k) = 0.0<br /> END DO<br /> END DO<br /> END DO<br /> END SUBROUTINE inicia_zeros<br /> <br /> SUBROUTINE initialize_V(V0,rho)<br /> ! Inicialização do potencial<br /> USE globais<br /> IMPLICIT NONE<br /> INTEGER :: i,j,k<br /> REAL,DIMENSION(L,L,L) :: V0,rho<br /> ! Inicializa o potencial <br /> DO k=1,L <br /> DO i=1,L<br /> DO j=1,L<br /> V0(j,i,k) = 0.0<br /> rho(j,i,k) = 0.0<br /> END DO<br /> END DO<br /> END DO<br /> rho(11,11,11) = 1.0/dx**3<br /> rho(11,11,10) = 1.0/dx**3<br /> rho(11,11,12) = 1.0/dx**3<br /> rho(11,10,11) = 1.0/dx**3<br /> rho(11,12,11) = 1.0/dx**3<br /> rho(10,11,11) = 1.0/dx**3<br /> rho(12,11,11) = 1.0/dx**3<br /> END SUBROUTINE initialize_V<br /> <br /> SUBROUTINE update_V_jacobi(Vn,Vn1,rho,deltaV)<br /> ! Calcula a atualização uma vez<br /> USE globais<br /> IMPLICIT NONE<br /> INTEGER :: i,j,k<br /> REAL :: deltaV<br /> REAL,DIMENSION(L,L,L) :: Vn,Vn1,rho<br /> deltaV = 0.0<br /> DO k=2,L-1 <br /> DO i=2,L-1<br /> DO j=2,L-1<br /> Vn1(j,i,k) = (Vn(j+1,i,k) + Vn(j-1,i,k) + &amp;<br /> Vn(j,i+1,k) + Vn(j,i-1,k) + Vn(j,i,k+1) + Vn(j,i,k-1))/6.0 &amp;<br /> + (rho(j,i,k)*dx**2)/6.0<br /> deltaV = deltaV + abs(Vn(j,i,k) - Vn1(j,i,k))<br /> END DO<br /> END DO<br /> END DO<br /> END SUBROUTINE update_V_jacobi <br /> <br /> SUBROUTINE laplace_jacobi(Vm,Vm1,rho,delta_v,cont)<br /> ! Implementa o método de Jacobi<br /> USE globais<br /> IMPLICIT NONE<br /> INTEGER :: cont<br /> REAL :: delta_v,tol<br /> REAL,DIMENSION(L,L,L) :: Vm,Vm1,rho<br /> tol = 1e-5*L*L*L <br /> delta_v = 100.0<br /> DO WHILE (delta_v .gt. tol)<br /> CALL update_V_jacobi(Vm,Vm1,rho,delta_v)<br /> CALL update_V_jacobi(Vm1,Vm,rho,delta_v)<br /> cont = cont+1<br /> END DO<br /> END SUBROUTINE laplace_jacobi<br /> <br /> SUBROUTINE update_V_gauss_siedel(Vn,rho,deltaV)<br /> ! Calcula a atualização uma vez<br /> USE globais<br /> IMPLICIT NONE<br /> INTEGER :: i,j,k<br /> REAL :: deltaV,v1,v2<br /> REAL,DIMENSION(L,L,L) :: Vn,rho<br /> deltaV = 0.0<br /> DO k=2,L-1<br /> DO i=2,L-1<br /> DO j=2,L-1<br /> v1 = Vn(j,i,k)<br /> Vn(j,i,k) = (Vn(j+1,i,k) + Vn(j-1,i,k) + &amp;<br /> Vn(j,i+1,k) + Vn(j,i-1,k) + Vn(j,i,k+1) + Vn(j,i,k-1))/6.0 &amp;<br /> + (rho(j,i,k)*dx**2)/6.0<br /> v2 = Vn(j,i,k)<br /> deltaV = deltaV + abs(v2 - v1)<br /> END DO<br /> END DO<br /> END DO<br /> END SUBROUTINE update_V_gauss_siedel<br /> <br /> SUBROUTINE laplace_gauss_siedel(Vm,rho,delta_v,cont)<br /> ! Implementa o método de Gauss-Siedel<br /> USE globais<br /> IMPLICIT NONE<br /> INTEGER :: cont<br /> REAL :: delta_v,tol<br /> REAL,DIMENSION(L,L,L) :: Vm,rho<br /> tol = 1e-5*L*L*L <br /> delta_v = 100.0<br /> DO WHILE (delta_v .gt. tol)<br /> CALL update_V_gauss_siedel(Vm,rho,delta_v)<br /> CALL update_V_gauss_siedel(Vm,rho,delta_v)<br /> cont = cont+1<br /> END DO<br /> END SUBROUTINE laplace_gauss_siedel<br /> <br /> SUBROUTINE electric_field(Vk,Ex,Ey,Ez)<br /> ! Calcula o campo elétrico a partir do potencial <br /> USE globais<br /> IMPLICIT NONE<br /> INTEGER :: i,j,k<br /> REAL,DIMENSION(L,L,L) :: Vk,Ex,Ey,Ez<br /> <br /> ! Bordas k=1 e k=L<br /> DO k=1,L,L-1<br /> DO i=1,L<br /> DO j=1,L<br /> ! Derivada para trás<br /> IF (i == L) THEN<br /> Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k))<br /> Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k))<br /> Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1))<br /> ! Derivada para trás<br /> ELSE IF (j == L) THEN<br /> Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k)) <br /> Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k)) <br /> Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1)) <br /> ! Derivada para trás<br /> ELSE IF (k == L) THEN<br /> Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k)) <br /> Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k)) <br /> Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1)) <br /> ! Derivada para frente<br /> ELSE<br /> Ex(j,i,k) = - (Vk(j,i+1,k) - Vk(j,i,k))<br /> Ey(j,i,k) = - (Vk(j+1,i,k) - Vk(j,i,k))<br /> Ez(j,i,k) = - (Vk(j,i,k+1) - Vk(j,i,k))<br /> END IF<br /> END DO<br /> END DO<br /> END DO<br /> ! Bordas j=1 e j=L<br /> DO j=1,L,L-1<br /> DO i=1,L<br /> DO k=1,L<br /> ! Derivada para trás<br /> IF (i == L) THEN<br /> Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k))<br /> Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k))<br /> Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1))<br /> ! Derivada para trás<br /> ELSE IF (j == L) THEN<br /> Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k)) <br /> Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k)) <br /> Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1)) <br /> ! Derivada para trás<br /> ELSE IF (k == L) THEN<br /> Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k)) <br /> Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k)) <br /> Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1)) <br /> <br /> ! Derivada para frente<br /> ELSE<br /> Ex(j,i,k) = - (Vk(j,i+1,k) - Vk(j,i,k))<br /> Ey(j,i,k) = - (Vk(j+1,i,k) - Vk(j,i,k))<br /> Ez(j,i,k) = - (Vk(j,i,k+1) - Vk(j,i,k))<br /> END IF<br /> END DO<br /> END DO<br /> END DO<br /> ! Bordas i=1 e i=L<br /> DO i=1,L,L-1<br /> DO k=1,L<br /> DO j=1,L<br /> ! Derivada para trás<br /> IF (i == L) THEN<br /> Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k))<br /> Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k))<br /> Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1))<br /> ! Derivada para trás<br /> ELSE IF (j == L) THEN<br /> Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k)) <br /> Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k)) <br /> Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1)) <br /> ! Derivada para trás<br /> ELSE IF (k == L) THEN<br /> Ex(j,i,k) = - (Vk(j,i,k) - Vk(j,i-1,k)) <br /> Ey(j,i,k) = - (Vk(j,i,k) - Vk(j-1,i,k)) <br /> Ez(j,i,k) = - (Vk(j,i,k) - Vk(j,i,k-1)) <br /> ! Derivada para frente<br /> ELSE<br /> Ex(j,i,k) = - (Vk(j,i+1,k) - Vk(j,i,k))<br /> Ey(j,i,k) = - (Vk(j+1,i,k) - Vk(j,i,k))<br /> Ez(j,i,k) = - (Vk(j,i,k+1) - Vk(j,i,k))<br /> END IF<br /> END DO<br /> END DO<br /> END DO<br /> <br /> ! Derivada centrada <br /> DO k=2,L-1<br /> DO i=2,L-1<br /> DO j=2,L-1 <br /> Ex(j,i,k) = - (Vk(j,i+1,k) - Vk(j,i-1,k))/2.0 <br /> Ey(j,i,k) = - (Vk(j+1,i,k) - Vk(j-1,i,k))/2.0 <br /> Ez(j,i,k) = - (Vk(j,i,k+1) - Vk(j,i,k-1))/2.0 <br /> END DO<br /> END DO<br /> END DO<br /> END SUBROUTINE electric_field<br /> <br /> SUBROUTINE dados(x1,y1,z1,Vk,Ex,Ey,Ez)<br /> ! Imprime os dados em arquivos<br /> USE globais<br /> IMPLICIT NONE<br /> INTEGER :: i,j,k,m<br /> REAL,DIMENSION(L,L,L) :: Vk,Ex,Ey,Ez,x1,y1,z1<br /> CHARACTER(28),DIMENSION(2) :: a<br /> CHARACTER(27),DIMENSION(2) :: b<br /> <br /> a(1) = 'potencial_poisson_jacob1.dat'<br /> a(2) = 'potencial_poisson_gauss1.dat'<br /> b(1) = 'electric_poisson_jacob1.dat'<br /> b(2) = 'electric_poisson_gauss1.dat'<br /> <br /> DO m=1,2<br /> OPEN(10,FILE=a(m),STATUS='unknown',ACTION='write')<br /> DO i=1,L<br /> DO j=1,L<br /> DO k=1,L<br /> WRITE(10,*) x1(j,i,k),y1(j,i,k),z1(j,i,k),Vk(j,i,k)<br /> END DO<br /> WRITE(10,*) ''<br /> END DO<br /> WRITE(10,*) ''<br /> END DO<br /> CLOSE(10)<br /> <br /> OPEN(20,FILE=b(m),STATUS='unknown',ACTION='write')<br /> DO i=1,L<br /> DO j=1,L<br /> DO k=1,L<br /> WRITE(20,*) x1(j,i,k),y1(j,i,k),z1(j,i,k), &amp;<br /> Ex(j,i,k),Ey(j,i,k),Ez(j,i,k)<br /> END DO<br /> WRITE(20,*) '' <br /> END DO<br /> WRITE(20,*) ''<br /> END DO<br /> CLOSE(20)<br /> END DO<br /> END SUBROUTINE dados</div>

Gustavobopsin