initial.f

二维热能方程的求解

!=======================================================================
      subroutine initial(nx,ny,idpsrc,level,id_hsrc,ax,bx,ay,by,at,x)
!=======================================================================
      implicit none
      integer nx,ny,idpsrc,level,id_hsrc
      real*8 ax,bx,ay,by,at
      real*8 x(nx,ny)

      integer ix,iy
      real*8 hx,hy,px,py,x1,y1,truesol

      if(level.ge.1) print'("INITIAL: @t=",f7.3)',at
      if(level.ge.2) then
         print*,"nx=",nx," ny=",ny
         print*,"idpsrc=",idpsrc," id_hsrc=",id_hsrc
         print*,"ax=",ax," bx=",bx
         print*,"ay=",ay," by=",by
      endif

      do iy=1,ny
      do ix=1,nx
         x(ix,iy)=0.0d0
      enddo
      enddo

      if(id_hsrc.ge.1) then
         if(level.ge.2) then
            print*,"INITIAL=0. id_hsrc=",id_hsrc
         endif
         return
      endif
      
      hx=(bx-ax)/dble(nx-1)
      hy=(by-ay)/dble(ny-1)
      x1=ax+(bx-ax)/5.
      y1=ay+(by-ay)/5.

      if(idpsrc.eq.1) then

         do iy=1,ny
            py=dble(iy-1)*hy+ay
         do ix=1,nx
            px=dble(ix-1)*hx+ax
            if(px.le.x1 .and. py.le.y1) then
               x(ix,iy)=1.0d0
            endif
         enddo
         enddo
         if(level.ge.2) then
            print*,"INITIAL=1.0 on the first of 5X5 subregions."
         endif

      else

         do iy=1,ny
            py=dble(iy-1)*hy+ay
         do ix=1,nx
            px=dble(ix-1)*hx+ax
            x(ix,iy)=truesol(at,px,py)
         enddo
         enddo
         if(level.ge.2) then
            print*,"INITIAL=truesol(t0,x,y)"
         endif

      endif

      return
      end