inputs.f

二维热能方程的求解

c=======================================================================
      real*8 function truesol(t,x,y)
c=======================================================================
      implicit none
      include 'commons.h'
      real*8 t,x,y

      if(idELL.eq.1)then
         truesol=dsin(freqx*x)*dsin(freqy*y)
         return
      endif

      if(ntrue.eq.0)then
         truesol=t+x*x+y*y
      else if(ntrue.eq.1)then
         truesol=dsin(freqt*t)+dsin(freqx*x)+dsin(freqy*y)
      else if(ntrue.eq.2)then
         truesol=dsin(freqt*t)*dsin(freqx*x)*dsin(freqy*y)
      else if(ntrue.eq.3)then
         truesol=dexp(-t)*dsin(freqx*x)*dsin(freqy*y)
      endif
      return
      end

c-----------------------------------------------------------------------
      real*8 function dtsol(t,x,y)
c-----------------------------------------------------------------------
      implicit none
      include 'commons.h'
      real*8 t,x,y

      if(idELL.eq.1)then
         dtsol=0.d0
         return
      endif

      if(ntrue.eq.0)then
         dtsol=1.d0
      else if(ntrue.eq.1)then
         dtsol=freqt*dcos(freqt*t)
      else if(ntrue.eq.2)then
         dtsol=freqt*dcos(freqt*t)*dsin(freqx*x)*dsin(freqy*y)
      else if(ntrue.eq.3)then
         dtsol=-dexp(-t)*dsin(freqx*x)*dsin(freqy*y)
      endif
      return
      end

c-----------------------------------------------------------------------
      real*8 function dxsol(t,x,y)
c-----------------------------------------------------------------------
      implicit none
      include 'commons.h'
      real*8 t,x,y

      if(idELL.eq.1)then
         dxsol=freqx*dcos(freqx*x)*dsin(freqy*y)
         return
      endif

      if(ntrue.eq.0)then
         dxsol=2.d0*x
      else if(ntrue.eq.1)then
         dxsol=freqx*dcos(freqx*x)
      else if(ntrue.eq.2)then
         dxsol=dsin(freqt*t)*freqx*dcos(freqx*x)*dsin(freqy*y)
      else if(ntrue.eq.3)then
         dxsol=dexp(-t)*freqx*dcos(freqx*x)*dsin(freqy*y)
      endif
      return
      end

c-----------------------------------------------------------------------
      real*8 function dysol(t,x,y)
c-----------------------------------------------------------------------
      implicit none
      include 'commons.h'
      real*8 t,x,y

      if(idELL.eq.1)then
         dysol=dsin(freqx*x)*freqy*dcos(freqy*y)
         return
      endif

      if(ntrue.eq.0)then
         dysol=2.d0*y
      else if(ntrue.eq.1)then
         dysol=freqy*dcos(freqy*y)
      else if(ntrue.eq.2)then
         dysol=dsin(freqt*t)*dsin(freqx*x)*freqy*dcos(freqy*y)
      else if(ntrue.eq.3)then
         dysol=dexp(-t)*dsin(freqx*x)*freqy*dcos(freqy*y)
      endif
      return
      end

c-----------------------------------------------------------------------
      real*8 function dxxsol(t,x,y)
c-----------------------------------------------------------------------
      implicit none
      include 'commons.h'
      real*8 t,x,y,truesol

      if(idELL.eq.1)then
         dxxsol=-freqx*freqx*dsin(freqx*x)*dsin(freqy*y)
         return
      endif

      if(ntrue.eq.0)then
         dxxsol=2.d0
      else if(ntrue.eq.1)then
         dxxsol=-freqx*freqx*dsin(freqx*x)
      else if(ntrue.eq.2)then
         dxxsol=-freqx*freqx*dsin(freqt*t)*dsin(freqx*x)*dsin(freqy*y)
      else if(ntrue.eq.3)then
         dxxsol=-freqx*freqx*truesol(t,x,y)
      endif
      return
      end

c-----------------------------------------------------------------------
      real*8 function dyysol(t,x,y)
c-----------------------------------------------------------------------
      implicit none
      include 'commons.h'
      real*8 t,x,y,truesol

      if(idELL.eq.1)then
         dyysol=-freqy*freqy*dsin(freqx*x)*dsin(freqy*y)
         return
      endif

      if(ntrue.eq.0)then
         dyysol=2.d0
      else if(ntrue.eq.1)then
         dyysol=-freqy*freqy*dsin(freqy*y)
      else if(ntrue.eq.2)then
         dyysol=-freqy*freqy*dsin(freqt*t)*dsin(freqx*x)*dsin(freqy*y)
      else if(ntrue.eq.3)then
         dyysol=-freqy*freqy*truesol(t,x,y)
      endif
      return
      end

c=======================================================================
      real*8 function b1(x,y)
c=======================================================================
      implicit none
      include 'commons.h'
      real*8 x,y
      if (ndiff.le.1) then
          b1=1.d0
      else if(ndiff.eq.2 .or. ndiff.eq.4)then
          b1=1.d0/(2.d0+dcos(3.d0*pi*x)*dcos(2.d0*pi*y))
      else if(ndiff.eq.3)then
          if(x.le.0.5001)then
            b1=1.d0+0.5*dsin(5.d0*pi*x)+y**3
          else
            b1=1.5d0/(1.d0+(x-0.5d0)**2)+y**3
          endif
      else
        stop 'Do definition for ndiff'
      endif
      return
      end

c-----------------------------------------------------------------------
      real*8 function dxb1(x,y)
c-----------------------------------------------------------------------
      implicit none
      include 'commons.h'
      real*8 x,y,b1
      if (ndiff.le.1) then
          dxb1=0.d0
      else if(ndiff.eq.2 .or. ndiff.eq.4)then
          dxb1=3.d0*pi*dsin(3.d0*pi*x)*dcos(2.d0*pi*y)*(b1(x,y)**2)
      else if(ndiff.eq.3)then
          if(x.le.0.5001)then
             dxb1=2.5d0*pi*dcos(5.d0*pi*x)
          else
             dxb1=-3.d0*(x-0.5d0)/(1.d0+(x-0.5d0)**2)**2
          endif
      endif
      return
      end

c=======================================================================
      real*8 function b2(x,y)
c=======================================================================
      implicit none
      include 'commons.h'
      real*8 x,y
      if (ndiff.le.1) then
          b2=1.d0
      else if(ndiff.eq.2)then
          b2=1.d0/(2.d0+dcos(3.d0*pi*x)*dcos(2.d0*pi*y))
      else if(ndiff.eq.3 .or. ndiff.eq.4)then
          if(x.le.0.5001)then
            b2=1.d0+0.5*dsin(5.d0*pi*x)+y**3
          else
            b2=1.5d0/(1.d0+(x-0.5d0)**2)+y**3
          endif
      endif
      return
      end

c-----------------------------------------------------------------------
      real*8 function dyb2(x,y)
c-----------------------------------------------------------------------
      implicit none
      include 'commons.h'
      real*8 x,y,b2
      if (ndiff.le.1) then
          dyb2=0.d0
      else if(ndiff.eq.2)then
          dyb2=2.d0*pi*dcos(3.d0*pi*x)*dsin(2.d0*pi*y)*(b2(x,y)**2)
      else if(ndiff.eq.3 .or. ndiff.eq.4)then
          if(x.le.0.5001)then
             dyb2=3.d0*y*y
          else
             dyb2=3.d0*y*y
          endif
      endif
      return
      end

c=======================================================================
      real*8 function ftnc(x,y)
c=======================================================================
      implicit none
      include 'commons.h'
      real*8 x,y
      if(nreact.eq.0)then
         ftnc=0.0d0
      else
         ftnc=0.0d0
      endif
      return
      end

c=======================================================================
      real*8 function ftncu(t,u)
c=======================================================================
      implicit none
      include 'commons.h'
      real*8 t,u

      if(nonlin.eq.0)then
         ftncu=0.0d0
      else
         ftncu=-(dexp(u)*u*u)/sqrt(0.1d0+t*t+2.d0*u*u)
      endif
      return
      end

c=======================================================================
      real*8 function beta(x,y,iside)
c=======================================================================
      implicit none
      include 'commons.h'
      integer iside
      real*8 x,y

      if(idELL.eq.1)then
         beta=100.d0
         return
      elseif(nbc.le.2) then
         beta=0.d0
         return
      endif

      goto(10,20,30,40),iside
 10   continue
         beta=0.d0
      return
 20   continue
         beta=0.d0
      return
 30   continue
         beta=0.d0
      return
 40   continue
         beta=0.d0
      return
      end

************************************************************************
***** RECOMMENDED NOT TO TOUCH BELOW ***********************************
************************************************************************

c=======================================================================
      real*8 function ftnf(t,x,y)
c=======================================================================
      implicit none
      include 'commons.h'
      real*8 t,x,y
      real*8 b1,b2,dxb1,dyb2,ftnc,ftncu,u
      real*8 truesol,dtsol,dxsol,dysol,dxxsol,dyysol

      u=truesol(t,x,y)

      ftnf=dtsol(t,x,y)+ftnc(x,y)*u-ftncu(t,u)
     &    -( dxb1(x,y)*dxsol(t,x,y)+b1(x,y)*dxxsol(t,x,y)
     &      +dyb2(x,y)*dysol(t,x,y)+b2(x,y)*dyysol(t,x,y) )
      return
      end

c=======================================================================
      real*8 function ftng(t,x,y,iside)
c=======================================================================
      implicit none
      include 'commons.h'
      integer iside
      real*8 t,x,y
      real*8 b1,b2,beta,truesol,dxsol,dysol
      real*8 bdry

      bdry=beta(x,y,iside)*truesol(t,x,y)

      goto(10,20,30,40),iside
 10   ftng=-b1(x,y)*dxsol(t,x,y)+bdry
      return
 20   ftng= b1(x,y)*dxsol(t,x,y)+bdry
      return
 30   ftng=-b2(x,y)*dysol(t,x,y)+bdry
      return
 40   ftng= b2(x,y)*dysol(t,x,y)+bdry
      return
      end