Ask Question

Name:
Title:
Your Question:

Answer Question

Name:
Your Answer:
User Submitted Source Code!


Description:
  df
Language: FORTRAN
Code:
! **********************************************************************
!  FICHERO: calor_euler_explicito.f
!  OTROS FICHEROS:
!  DESCRICION:  
!    resolvemos la ecuacion del calor utilizando el metodo
!    Euler Explicito para discretizacion temporal y
!    Diferencias Finitas para discretizacion en espacio
!  AUTOR: Pedro Galan del Sastre
!  ULTIMA REVISION: 30/04/14
! **********************************************************************

      program ecuacion_calor

      integer   i, j, n, npt
      integer   nx, ny
      real*8    a, b, c, d
      real*8    hx, hy
      real*8    nu, dt
      real*8    x(*), y(*), u(*), un(*)
      real*8    cc_i(*), cc_d(*), cc_ab(*), cc_ar(*)
      pointer   (ptr_x, x)
      pointer   (ptr_y, y)
      pointer   (ptr_u, u)
      pointer   (ptr_un, un)
      pointer   (ptr_cc_i, cc_i)
      pointer   (ptr_cc_d, cc_d)
      pointer   (ptr_cc_ab, cc_ab)
      pointer   (ptr_cc_ar, cc_ar)
      integer   ind


! definimos el dominio
      a = 0d0
      b = 1d0
      c = 0d0
      d = 1d0

! variables del problema
      hx = 0.05d0
      hy = 0.05d0
      nx = 1 + 1d0/hx
      ny = 1 + 1d0/hy
      nu = 0.001d0
      dt = 0.5d0
      npt = 300d0/dt

! reservamos memoria
      ptr_x = malloc(nx*ny*8)
      ptr_y = malloc(nx*ny*8)
      ptr_u = malloc(nx*ny*8)
      ptr_un = malloc(nx*ny*8)
      ptr_cc_i = malloc(ny*8)
      ptr_cc_d = malloc(ny*8)
      ptr_cc_ab = malloc(nx*8)
      ptr_cc_ar = malloc(nx*8)

! definimos coordenadas de los nodos
      do 10 i=1,nx
        do 20 j=1,ny
           x(ind(i,j,ny)) = a + hx*i
           y(ind(i,j,ny)) = c + hy*i
 20     continue
 10   continue

! valor inicial y condiciones de contorno
      do 30 i=1,nx
        do 40 j=1,ny
           un(ind(i,j,ny)) = 30
 40     continue
 30   continue

      do 50 i=1,nx
        cc_ab(i) = 0d0
        cc_ar(i) = 50d0
 50   continue

      do 60 i=1,ny
        cc_i(i) = 0d0
        cc_d(i) = 0d0
 60   continue


! bucle con los pasos de tiempo
      open (unit=11, file='u.dat', action='write', 
     &status='replace', access="stream", form='unformatted')
      write(11) (un(i),i=1,nx*ny)

      do 100 n=1,npt
        call euler_explicito_unpaso(u, un, nu, dt,
     &                  cc_i, cc_d, cc_ab, cc_ar, hx, hy, nx, ny)

        write(11) (u(i),i=1,nx*ny)

        do 110 i=1,nx
           do 120 j=1,ny
                un(ind(i,j,ny)) = u(ind(i,j,ny))
 120       continue
 110            continue

 100  continue

! liberamos memoria
      call free(ptr_x);
      call free(ptr_y);
      call free(ptr_u);
      call free(ptr_un);
      call free(ptr_cc_i);
      call free(ptr_cc_d);
      call free(ptr_cc_ab);
      call free(ptr_cc_ar);

      end


! **************************************************************************
! funcion euler_explicito_unpaso
! DESCRICION:
!   Calculamos un paso de tiempo del esquema Euler Explicito
!   utilizando Diferencias Finitas para discretizacion en tiempo
! **************************************************************************

      subroutine euler_explicito_unpaso(u, un, nu, dt,
     &                  cc_i, cc_d, cc_ab, cc_ar, hx, hy, nx, ny)
      integer   i, j
      integer   nx, ny
      real*8    hx, hy
      real*8    nu, dt
      real*8    u(*), un(*)
      real*8    cc_i(*), cc_d(*), cc_ab(*), cc_ar(*)
      integer   ind

      do 210 i=2,nx-1
        do 220 j=2,ny-1
           u(ind(i,j,ny)) = un(ind(i,j,ny)) + nu*dt*(un(ind(i+1,j,ny)) +
     &                  un(ind(i-1,j,ny)) - 4d0*un(ind(i,j,ny)) +
     &                  un(ind(i,j+1,ny)) + un(ind(i,j-1,ny)))/(hx*hy)
 220    continue
 210  continue

      do 310 i=1,nx
        u(ind(i,1,ny)) = cc_ab(i)
        u(ind(i,ny,ny)) = cc_ar(i)
 310  continue

      do 320 i=1,ny
        u(ind(1,i,ny)) = cc_i(i)
        u(ind(nx,i,ny)) = cc_d(i)
 320  continue

      end

! **************************************************************************
! funcion ind
! DESCRICION:
!   Calcula el nodo k con  indices (i,j)
! **************************************************************************

      integer function ind(i,j,ny)

      ind = (i-1)*ny + (j-1) + 1

      end
          
          
Comments: