module utils implicit none contains ! Drills a hole of radius r in the array holes centered in (x,y) ! and sets a thermostat at temperature t in this hole subroutine drill_hole(holes, thermostats, nx, ny, x, y, delta, r, t) implicit none real(8), dimension(:,:), allocatable :: holes, thermostats real(8) :: delta real :: x,y,r,t integer :: nx, ny, i, j do i = 0, nx - 1 do j = 0, ny - 1 if ((float(i) * delta - x)**2.0 + (float(j) * delta - y)**2.0 <= r**2.0) then holes(i,j) = 0 thermostats(i,j) = t endif enddo enddo end end ! Programme principal program main use utils implicit none real(8) :: eps, delta, x, y integer :: niter_max, p, nx, ny, i, j real(8), dimension(:,:), allocatable :: Tp, Tn, Q, GradX, GradY, holes niter_max = 10000 eps = 1e-5 nx = 100 ny = 100 delta = 0.01 allocate(Tp(0:nx-1,0:ny-1)) allocate(Tn(0:nx-1,0:ny-1)) allocate(Q(0:nx-1,0:ny-1)) allocate(GradX(0:nx-1,0:ny-1)) allocate(GradY(0:nx-1,0:ny-1)) allocate(holes(0:nx-1,0:ny-1)) ! Initialize Q do i = 0, nx-1 do j = 0, ny-1 Q(i, j) = 0 enddo enddo ! Initialize holes do i = 0, nx-1 do j = 0, ny-1 holes(i,j) = 1 enddo enddo ! First mickey call drill_hole(holes, Q, nx, ny, 0.30, 0.50, delta, 0.10, 0.5) call drill_hole(holes, Q, nx, ny, 0.20, 0.60, delta, 0.02, 1.0) call drill_hole(holes, Q, nx, ny, 0.40, 0.60, delta, 0.02, 0.0) ! Second mickey call drill_hole(holes, Q, nx, ny, 0.70, 0.50, delta, 0.10, 0.5) call drill_hole(holes, Q, nx, ny, 0.60, 0.60, delta, 0.02, 1.0) call drill_hole(holes, Q, nx, ny, 0.80, 0.60, delta, 0.02, 0.0) ! Initialize side conditions Tp = Q ! Gauss-Seidel loop do p = 1, niter_max ! Iterate do i = 1, nx-1 do j = 1, ny-1 if (holes(i,j) == 1) then Tn(i,j) = (delta * delta * Q(i,j) + Tp(i+1,j) + Tp(i,j+1) + Tn(i-1,j) + Tn(i,j-1)) / 4 else Tn(i,j) = Q(i,j) endif enddo enddo ! Exit if the algorithm converged if (maxval(abs(Tp-Tn)) <= eps) exit Tp = Tn enddo Tn = Tp ! Gradient computation do i = 0, nx-1 do j = 0, ny-1 if (i == 0) then GradX(i,j) = (Tn(i+1,j) - Tn(i ,j)) / delta elseif (i == nx-1) then GradX(i,j) = (Tn(i,j ) - Tn(i-1,j)) / delta else GradX(i,j) = (Tn(i+1,j) - Tn(i-1,j)) / (2 * delta) endif if (j == 0) then GradY(i,j) = (Tn(i,j+1) - Tn(i,j )) / delta elseif (j == ny-1) then GradY(i,j) = (Tn(i,j ) - Tn(i,j-1)) / delta else GradY(i,j) = (Tn(i,j+1) - Tn(i,j-1)) / (2 * delta) endif enddo enddo write(*,*) 'Iterations : ', p ! Output open(1, file='culasse.res') do i = 0, nx-1 do j = 0, ny-1 write(1,*) (float(i) * delta), ' ', (float(j) * delta), ' ', Tn(i,j), ' ', GradX(i,j), ' ', GradY(i,j) enddo enddo close(1) deallocate(Tp) deallocate(Tn) deallocate(Q) deallocate(GradX) deallocate(GradY) deallocate(holes) end