module resultados ! objetivo: calcular solução numérica e apresentar gráficos dos resultados !------------------------------------------------- use coeficientes use solvers use portlib implicit none contains !------------------------------------------------- subroutine norma_L1 ! subroutine norma_L1 (nx,ny,nxy,a,b,x,norma) ! integer,intent(in) :: nx, ny, nxy ! real*8,intent(in),dimension(nxy,5) :: a ! coeficientes ! real*8,intent(in),dimension(nxy) :: b ! fonte ! real*8,intent(in),dimension(nxy) :: x ! solução ! real*8,intent(inout) :: norma norma = 0.0d0 do i = 2, nx-1 do j = 2, ny-1 np = (j-1)*nx + i norma = norma + dabs(T(np)*ap - as*T(np-nx) & - aw*T(np-1) - ae*T(np+1) - an*T(np+nx) ) ! norma = norma + dabs(x(np)*a(np,3) + a(np,1)*x(np-nx) & ! + a(np,2)*x(np-1) + a(np,4)*x(np+1) & ! + a(np,5)*x(np+nx) - b(np)) end do end do end subroutine norma_L1 !------------------------------------------------- subroutine solucao_numerica_monogrid ! real*8 :: norma_1 ! norma da primeira iteração tcpu = timef() ! alocação de memória allocate ( t(nxyg(1)) ) ! allocate ( t(nxyg(1)), a(nxyg(1),5), b(nxyg(1)) ) ! if ( solver == 5 ) allocate ( dl(nxyg(1),4), du(nxyg(1),3) ) ! t = 0.0d0 nx = nxg(1) ny = nyg(1) nxy = nx*ny ! nmed = (ny/2)*nx + nx nmed = (ny/2)*nx + nx/2 + 1 ! dx = 5.0d-1 / (nx-1) dx = 1.0d+0 / (nx-1) dy = 1.0d+0 / (ny-1) ! call inicializa_contornos ! call inicializa_analitico if ( tipo_ci == 0) t = 0.0d0 if ( tipo_ci == 1) call inicializa_analitico call calcula_T_contornos call calcula_coeficientes ! call norma_L1 (nx,ny,nxy,a,b,t,norma_i) ! if ( solver == 5 ) call lu2d5(a,dl,du,nxy,nx,ny) do it = 1, itmax ! início do ciclo iterativo ! resolve T da malha atual ! if ( solver == 1 ) call jacobi (nx,ny,nxy,nitm,a,b,t) ! if ( solver == 2 ) call gauss_seidel (nx,ny,nxy,nitm,a,b,t) ! call gauss_seidel (nx,ny,nxy,nitm,a,b,t) ! if ( solver == 3 ) call sor (nx,ny,nxy,nitm,w,a,b,t) ! if ( solver == 4 ) call adi (nx,ny,nxy,nitm,a,b,t) ! if ( solver == 5 ) call fb2d5(a,dl,du,ny,nx,nxy,t,tol,nitm,b) call resolve_T ! call norma_L1 (nx,ny,nxy,a,b,t,norma) call norma_L1 if ( it == 1 ) norma_i = norma ! calcula e escreve variáveis de interesse if ( it==1 .or. it==itmax .or. mod(it,freq)==0) then call globais write(8,2) it, t(nmed), tmed 2 format(i8,2x,2(1pe24.15)) write(10,4) it, norma/norma_i 4 format(2x,i8,1pe15.6) end if if ( norma/norma_i <= tolt ) exit end do ! fim do ciclo iterativo tcpu = timef() write(10,6) norma_i 6 format(/,1pe15.6,' = norma inicial') if ( escreve == 1 ) then write(10,5) 5 format(//,30x,'TEMPERATURAS') call print1(t,ny,nx,nxy) end if end subroutine solucao_numerica_monogrid !------------------------------------------------- subroutine solucao_numerica_multigrid ! real*8 :: norma_1 ! norma da primeira iteração integer :: l, ic, jc, npac, i1, j1, npa1, npga, npgf, npf ! auxiliares tcpu = timef() ! allocate ( t(nxyg(nm)), r(npt), dt(npt) ) ! allocate ( t(nxyg(nm)), r(npt), dt(npt-nxyg(nm)) ) allocate ( t(nxyg(nm)), r(npt-nxyg(nm)), dt(npt-nxyg(nm)) ) ! t = 0.0d0 r = 0.0d0 dt = 0.0d0 m = nm ! solução inicial da malha mais fina npa = 0 do l = 1, m-1 npa = npa + nxyg(l) end do ! allocate ( a(nxyg(m),5), b(nxyg(m)) ) ! if ( solver == 5 ) allocate ( dl(nxyg(m),4), du(nxyg(m),3) ) ! inicialização de parâmetros da malha atual nx = nxg(m) ny = nyg(m) nxy = nx*ny ! dx = 5.0d-1 / (nx-1) dx = 1.0d+0 / (nx-1) dy = 1.0d+0 / (ny-1) if ( tipo_ci == 0) t = 0.0d0 if ( tipo_ci == 1) call inicializa_analitico if ( itmax == 0 ) go to 7 call calcula_T_contornos call calcula_coeficientes ! call norma_L1 (nx,ny,nxy,a,b,t,norma_i) ! if ( solver == 5 ) call lu2d5(a,dl,du,nxy,nx,ny) ! resolve T da malha atual ! if ( solver == 1 ) call jacobi (nx,ny,nxy,nitm,a,b,t) ! if ( solver == 2 ) call gauss_seidel (nx,ny,nxy,nitm,a,b,t) ! call gauss_seidel (nx,ny,nxy,nitm,a,b,t) ! if ( solver == 3 ) call sor (nx,ny,nxy,nitm,w,a,b,t) ! if ( solver == 4 ) call adi (nx,ny,nxy,nitm,a,b,t) ! if ( solver == 5 ) call fb2d5(a,dl,du,ny,nx,nxy,t,tol,nitm,b) call resolve_T call norma_L1 norma_i = norma ! call norma_L1 (nx,ny,nxy,a,b,t,norma_i) ! write(10,4) 1, norma/norma_i write(10,4) 1, 1.0d0 call globais ! nmed = (ny/2)*nx + nx nmed = (ny/2)*nx + nx/2 + 1 write(8,2) 1, t(nmed), tmed ! cálcula resíduo da malha mais fina ! do i = 2, nx-1 ! do j = 2, ny-1 ! np = (j-1)*nx + i ! npg = npa + (j-1)*nxg(m) + i ! r(npg) = - ap*T(np) + as*T(np-nx) & ! + aw*T(np-1) + ae*T(np+1) & ! + an*T(np+nx) ! r(npg) = - a(np,3)*t(np) - a(np,1)*t(np-nx) & ! - a(np,2)*t(np-1) - a(np,4)*t(np+1) & ! - a(np,5)*t(np+nx) + b(np) ! end do ! end do ! cálcula resíduo da malha mais fina e o injeta na malha grossa adjacente nx = nxg(nm-1) ny = nyg(nm-1) nxy = nx*ny npa = 0 do l = 1, nm-2 npa = npa + nxyg(l) end do do i = 2, nx-1 i1 = 2*(i-1) + 1 ! i malha mais fina do j = 2, ny-1 j1 = 2*(j-1) + 1 ! j malha mais fina npf = (j1-1)*nxg(nm) + i1 ! malha mais fina npg = npa + (j-1)*nxg(nm-1) + i r(npg) = - ap*T(npf) + as*T(npf-nxg(nm)) & + aw*T(npf-1) + ae*T(npf+1) & + an*T(npf+nxg(nm)) end do end do ! deallocate ( a, b ) ! if ( solver == 5 ) deallocate ( dl, du ) do it = 2, itmax ! início do ciclo iterativo ! ciclo Multigrid de restrição (malha fina -> grossa) do m = nm-1, 1, -1 npa = 0 do l = 1, m-1 npa = npa + nxyg(l) end do ! allocate ( dtl(nxyg(m)) ) ! allocate ( dtl(nxyg(m)), a(nxyg(m),5), b(nxyg(m)) ) ! if ( solver == 5 ) allocate ( dl(nxyg(m),4), du(nxyg(m),3) ) nx = nxg(m) ny = nyg(m) nxy = nx*ny ! dx = 5.0d-1 / (nx-1) dx = 1.0d+0 / (nx-1) dy = 1.0d+0 / (ny-1) call calcula_coeficientes ! b = 0.0d0 if ( m /= (nm-1) ) then ! Restrição (pontos internos) npa1 = npa + nxyg(m) do i = 2, nx-1 i1 = 2*(i-1) + 1 ! i malha fina do j = 2, ny-1 ! np = (j-1)*nx + i j1 = 2*(j-1) + 1 ! j malha fina npgf = npa1 + (j1-1)*nxg(m+1) + i1 npg = npa + (j-1)*nxg(m) + i r(npg) = r(npgf) ! b(np) = r(npg) end do end do end if ! dtl = 0.0d0 do i = 1, nx do j = 1, ny npg = npa + (j-1)*nxg(m) + i dt(npg) = 0.0d0 end do end do ! if ( solver == 5 ) call lu2d5(a,dl,du,nxy,nx,ny) ! resolve dtl da malha atual ! if ( solver == 1 ) call jacobi (nx,ny,nxy,nitm,a,b,dtl) ! if ( solver == 2 ) call gauss_seidel (nx,ny,nxy,nitm,a,b,dtl) ! call gauss_seidel (nx,ny,nxy,nitm,a,b,dtl) ! if ( solver == 3 ) call sor (nx,ny,nxy,nitm,w,a,b,dtl) ! if ( solver == 4 ) call adi (nx,ny,nxy,nitm,a,b,dtl) ! if ( solver == 5 ) call fb2d5(a,dl,du,ny,nx,nxy,dtl,tol,nitm,b) ! call resolve_dtl call resolve_dt ! atualiza dt global com solução da malha atual ! do i = 1, nx ! do j = 1, ny ! np = (j-1)*nx + i ! npg = npa + (j-1)*nxg(m) + i ! dt(npg) = dtl(np) ! end do ! end do do i = 2, nx-1 do j = 2, ny-1 ! np = (j-1)*nx + i npg = npa + (j-1)*nxg(m) + i r(npg) = - ap*dt(npg) + as*dt(npg-nx) & + aw*dt(npg-1) + ae*dt(npg+1) & + an*dt(npg+nx) + r(npg) ! r(npg) = - ap*dtl(np) + as*dtl(np-nx) & ! + aw*dtl(np-1) + ae*dtl(np+1) & ! + an*dtl(np+nx) + r(npg) ! r(npg) = - a(np,3)*dtl(np) - a(np,1)*dtl(np-nx) & ! - a(np,2)*dtl(np-1) - a(np,4)*dtl(np+1) & ! - a(np,5)*dtl(np+nx) + b(np) end do end do ! deallocate ( dtl ) ! deallocate ( dtl, a, b ) ! if ( solver == 5 ) deallocate ( dl, du ) end do ! fim do ciclo de restrição ! ciclo Multigrid de prolongação (malha grossa -> fina) do m = 2, nm npa = 0 do l = 1, m-1 npa = npa + nxyg(l) end do ! if ( solver == 5 ) allocate ( dl(nxyg(m),4), du(nxyg(m),3) ) nx = nxg(m) ny = nyg(m) nxy = nx*ny ! dx = 5.0d-1 / (nx-1) dx = 1.0d+0 / (nx-1) dy = 1.0d+0 / (ny-1) ! if ( m == nm ) then ! do i = 1, nx ! do j = 1, ny ! np = (j-1)*nx + i ! npg = npa + (j-1)*nxg(m) + i ! dt(npg) = 0.0d0 ! end do ! end do ! end if ! Prolongação (pontos internos) ! npac = npa - nxyg(m-1) ! do i = 2, nx-1 ! ic = (i-1)/2 + 1 ! i malha grossa (anterior) ! do j = 2, ny-1 ! np = (j-1)*nx + i ! jc = (j-1)/2 + 1 ! j malha grossa (anterior) ! npga = npac + (jc-1)*nxg(m-1) + ic ! npg = npa + (j-1)*nxg(m) + i ! if ( mod(i,2) == 0 ) then ! i par ! if ( mod(j,2) == 0 ) then ! j par ! dt(npg) = 2.5d-1*(dt(npga)+dt(npga+1)+dt(npga+nxg(m-1)) & ! + dt(npga+1+nxg(m-1))) + dt(npg) ! else ! j ímpar ! dt(npg) = 5.0d-1*(dt(npga)+dt(npga+1)) + dt(npg) ! end if ! else ! i ímpar ! if ( mod(j,2) == 0 ) then ! j par ! dt(npg) = 5.0d-1*(dt(npga)+dt(npga+nxg(m-1))) + dt(npg) ! else ! j ímpar ! dt(npg) = dt(npga) + dt(npg) ! end if ! end if ! end do ! end do if ( m == nm ) then ! allocate ( a(nxyg(m),5), b(nxyg(m)) ) ! Prolongação (pontos internos) npac = npa - nxyg(m-1) do i = 2, nx-1 ic = (i-1)/2 + 1 ! i malha grossa (anterior) do j = 2, ny-1 np = (j-1)*nx + i jc = (j-1)/2 + 1 ! j malha grossa (anterior) npga = npac + (jc-1)*nxg(m-1) + ic npg = npa + (j-1)*nxg(m) + i if ( mod(i,2) == 0 ) then ! i par if ( mod(j,2) == 0 ) then ! j par T(np) = 2.5d-1*(dt(npga)+dt(npga+1)+dt(npga+nxg(m-1)) & + dt(npga+1+nxg(m-1))) + T(np) else ! j ímpar T(np) = 5.0d-1*(dt(npga)+dt(npga+1)) + T(np) end if else ! i ímpar if ( mod(j,2) == 0 ) then ! j par T(np) = 5.0d-1*(dt(npga)+dt(npga+nxg(m-1))) + T(np) else ! j ímpar T(np) = dt(npga) + T(np) end if end if end do end do ! do i = 1, nx ! do j = 1, ny ! np = (j-1)*nx + i ! npg = npa + (j-1)*nxg(m) + i ! t(np) = t(np) + dt(npg) ! end do ! end do ! call inicializa_contornos call calcula_coeficientes ! if ( solver == 5 ) call lu2d5(a,dl,du,nxy,nx,ny) ! resolve T da malha mais fina ! if ( solver == 1 ) call jacobi (nx,ny,nxy,nitm,a,b,t) ! if ( solver == 2 ) call gauss_seidel (nx,ny,nxy,nitm,a,b,t) ! call gauss_seidel (nx,ny,nxy,nitm,a,b,t) ! if ( solver == 3 ) call sor (nx,ny,nxy,nitm,w,a,b,t) ! if ( solver == 4 ) call adi (nx,ny,nxy,nitm,a,b,t) ! if ( solver == 5 ) call fb2d5(a,dl,du,ny,nx,nxy,t,tol,nitm,b) call resolve_T call norma_L1 ! call norma_L1 (nx,ny,nxy,a,b,t,norma) ! escreve variáveis de interesse da malha mais fina if ( it==itmax .or. mod(it,freq)==0) then call globais ! nmed = (ny/2)*nx + nx nmed = (ny/2)*nx + nx/2 + 1 write(8,2) it, t(nmed), tmed write(10,4) it, norma/norma_i 4 format(2x,i8,1pe15.6) end if if ( norma/norma_i <= tolt ) then ! tcpu = timef() ! write(10,5) ! 5 format(//,30x,'TEMPERATURAS') ! call print1(t,ny,nx,nxy) ! deallocate ( a, b ) ! if ( solver == 5 ) deallocate ( dl, du ) goto 7 end if ! cálcula resíduo da malha mais fina ! do i = 2, nx-1 ! do j = 2, ny-1 ! np = (j-1)*nx + i ! npg = npa + (j-1)*nxg(m) + i ! r(npg) = - ap*T(np) + as*T(np-nx) & ! + aw*T(np-1) + ae*T(np+1) & ! + an*T(np+nx) ! r(npg) = - a(np,3)*t(np) - a(np,1)*t(np-nx) & ! - a(np,2)*t(np-1) - a(np,4)*t(np+1) & ! - a(np,5)*t(np+nx) + b(np) ! end do ! end do ! cálcula resíduo da malha mais fina e o injeta na malha grossa adjacente nx = nxg(nm-1) ny = nyg(nm-1) nxy = nx*ny npa = 0 do l = 1, nm-2 npa = npa + nxyg(l) end do do i = 2, nx-1 i1 = 2*(i-1) + 1 ! i malha mais fina do j = 2, ny-1 j1 = 2*(j-1) + 1 ! j malha mais fina npf = (j1-1)*nxg(nm) + i1 ! malha mais fina npg = npa + (j-1)*nxg(nm-1) + i r(npg) = - ap*T(npf) + as*T(npf-nxg(nm)) & + aw*T(npf-1) + ae*T(npf+1) & + an*T(npf+nxg(nm)) end do end do ! deallocate ( a, b ) else ! Prolongação (pontos internos) npac = npa - nxyg(m-1) do i = 2, nx-1 ic = (i-1)/2 + 1 ! i malha grossa (anterior) do j = 2, ny-1 np = (j-1)*nx + i jc = (j-1)/2 + 1 ! j malha grossa (anterior) npga = npac + (jc-1)*nxg(m-1) + ic npg = npa + (j-1)*nxg(m) + i if ( mod(i,2) == 0 ) then ! i par if ( mod(j,2) == 0 ) then ! j par dt(npg) = 2.5d-1*(dt(npga)+dt(npga+1)+dt(npga+nxg(m-1)) & + dt(npga+1+nxg(m-1))) + dt(npg) else ! j ímpar dt(npg) = 5.0d-1*(dt(npga)+dt(npga+1)) + dt(npg) end if else ! i ímpar if ( mod(j,2) == 0 ) then ! j par dt(npg) = 5.0d-1*(dt(npga)+dt(npga+nxg(m-1))) + dt(npg) else ! j ímpar dt(npg) = dt(npga) + dt(npg) end if end if end do end do ! allocate ( dtl(nxyg(m)) ) ! allocate ( dtl(nxyg(m)), a(nxyg(m),5), b(nxyg(m)) ) call calcula_coeficientes ! b = 0.0d0 ! reinicialização dos resíduos e dt ! do i = 2, nx-1 ! do j = 2, ny-1 ! np = (j-1)*nx + i ! npg = npa + (j-1)*nxg(m) + i ! b(np) = r(npg) ! end do ! end do ! do i = 1, nx ! do j = 1, ny ! np = (j-1)*nx + i ! npg = npa + (j-1)*nxg(m) + i ! dtl(np) = dt(npg) ! end do ! end do ! if ( solver == 5 ) call lu2d5(a,dl,du,nxy,nx,ny) ! resolve dtl da malha atual ! if ( solver == 1 ) call jacobi (nx,ny,nxy,nitm,a,b,dtl) ! if ( solver == 2 ) call gauss_seidel (nx,ny,nxy,nitm,a,b,dtl) ! call gauss_seidel (nx,ny,nxy,nitm,a,b,dtl) ! if ( solver == 3 ) call sor (nx,ny,nxy,nitm,w,a,b,dtl) ! if ( solver == 4 ) call adi (nx,ny,nxy,nitm,a,b,dtl) ! if ( solver == 5 ) call fb2d5(a,dl,du,ny,nx,nxy,dtl,tol,nitm,b) ! call resolve_dtl call resolve_dt ! atualiza dt global com solução da malha atual ! do i = 1, nx ! do j = 1, ny ! np = (j-1)*nx + i ! npg = npa + (j-1)*nxg(m) + i ! dt(npg) = dtl(np) ! end do ! end do ! deallocate ( dtl ) ! deallocate ( dtl, a, b ) end if ! if ( solver == 5 ) deallocate ( dl, du ) end do ! fim do ciclo de prolongação end do ! fim do ciclo iterativo 7 continue tcpu = timef() write(10,6) norma_i 6 format(/,1pe15.6,' = norma inicial') if ( itmax == 0 ) then call globais ! nmed = (ny/2)*nx + nx nmed = (ny/2)*nx + nx/2 + 1 end if if ( escreve == 1 ) then write(10,5) 5 format(//,30x,'TEMPERATURAS') call print1(t,ny,nx,nxy) end if 2 format(i8,2x,2(1pe24.15)) end subroutine solucao_numerica_multigrid !------------------------------------------------- subroutine globais tmed = 0.0d0 do i = 1, nx-1 do j = 1, ny-1 np = (j-1)*nx + i tmed = tmed + 2.5d-1*(t(np)+t(np+1)+t(np+nx)+t(np+nx+1)) end do end do tmed = tmed / ((nx-1)*(ny-1)) end subroutine globais !------------------------------------------------- subroutine print1 ( b, nj, ni, nij ) integer,intent(in) :: nij, ni, nj real*8,intent(in),dimension(nij) :: b ! incógnita a imprimir integer :: k, ll, ipart, npi, nps, l k=0 ll=8 ipart = 1 2 write(10,11)ipart do 1 j=1,nj np=(j-1)*ni + k npi=np+1 nps=np+ll if(ni.lt.8) nps=np+ni write(10,10)(b(l),l=npi,nps) 1 continue if(ni.lt.8) goto 3 if(nps.eq.nij) goto 3 k=k+8 ll=ni-k if(ll.gt.8)ll=8 ipart = ipart + 1 goto 2 10 format(8(1pe15.6)) 11 format(/,2x,'parte ',i3) 3 continue end subroutine print1 !------------------------------------------------- subroutine rex_numerico (arquivo,t,tmed) character*20,intent(in) :: arquivo real*8,intent(in),dimension(nxyg(nm)) :: t ! temperatura real*8,intent(in) :: tmed ! T média open(7,file=arquivo) ! número total de variáveis = 2 if ( delta == 1 ) then write(7,1) dx, dia, 2 1 format(1pe30.19,' = dx',10x,'Dia = ',a12,//, & i6,' = número total de variáveis',/) else write(7,2) dy, dia, 2 2 format(1pe30.19,' = dy',10x,'Dia = ',a12,//, & i6,' = número total de variáveis',/) end if ! parâmetro 1 (local) = temperatura central write(7,3) title, t(nmed) 3 format(3x,"'Temperatura Central, ",a50,"'",/, & 1pe30.19,' = solução numérica',/) ! parâmetro 2 (global) = temperatura média do campo write(7,4) title, tmed 4 format(3x,"'T Média, ",a50,"'",/, & 1pe30.19,' = solução numérica',/) close(7) end subroutine rex_numerico ! ----------------------------------------------- end module resultados