heat_trans.c

二维热能方程的求解

#include <stdio.h>
#include <math.h>
#include <time.h>
typedef  double ATYPE;
typedef  double VTYPE;
typedef  float DTYPE;

void heat_trans(nt,nx,ny,itmax,idELL,idsol,id_nl,idpcond,level,
            maxADI,maxPCG,maxGS,idSGS,
         ndiff,nonlin,nreact,nbc,ntrue,idpsrc,
            n1,n2,interpol,id_hsrc,idsymm,
         at,bt,ax,bx,ay,by,tol,wavet,wavex,wavey,tsbgn,tsend,eta,
         K4,A3,A5,A,x,b,src,wksp,ws,HEATSRC, ierr)

int   nt,nx,ny,itmax,idELL,idsol,id_nl,idpcond,level,maxADI,maxPCG,maxGS,idSGS;
int   ndiff,nonlin,nreact,nbc,ntrue,idpsrc,n1,n2,interpol,id_hsrc,idsymm;
VTYPE *at,*bt,*ax,*bx,*ay,*by,*tol,*wavet,*wavex,*wavey,*tsbgn,*tsend,*eta;
ATYPE A3[][1];            /* Is has dimension ny*nx*(2*nx+1) or */
ATYPE K4[];
ATYPE A5[];               /* Is has dimension ny*nx*(2*nx+1) or */
ATYPE A[][nx*ny][3];      /* First component has dimension 4. */
ATYPE x[][nx*ny];         /* First component has dimension 2. */
ATYPE b[];                /* This has dimension nx*ny. */
ATYPE src[][nx*ny];
ATYPE wksp[],ws[];
DTYPE HEATSRC[];
int   *ierr;

{
   int i,j,imode,ida1,ida2,ixlast,level0;
   ATYPE coefI,coefA1,coefA2;
   float utime0,stime0,utime1,stime1;

   if(level>=1)
       printf("HEAT_TRANS: nt=%d nx=%d ny=%d id_nl=%d\n",nt,nx,ny,id_nl);
   if(level>=1)
       printf(" ntrue=%d ndiff=%d nonlin=%d nreact=%d nbc=%d tol=%.3g eta=%g\n",
                     ntrue,ndiff,nonlin,nreact,nbc,*tol,*eta);
   if(level>=1) printf(" wavet=%g wavex=%g wavey=%g\n",*wavet,*wavex,*wavey);

   if(level>=3){
      printf(" at=%g bt=%g ax=%g bx=%g ay=%g by=%g\n",*at,*bt,*ax,*bx,*ay,*by);
      printf(" itmax=%d tol=%g idsol=%d idpcond=%d\n",itmax,*tol,idsol,idpcond);
      printf(" n1=%d n2=%d",n1,n2);
      printf(" interpol=%d id_hsrc=%d idsymm=%d\n",interpol,id_hsrc,idsymm);
      printf(" tsbgn=%g tsend=%g\n",*tsbgn,*tsend);
      }

   if(level>=1){
      #if defined (_IBMR2) || defined (AIX) || defined (hpux)
         etimef77(&utime0,&stime0);
      #else
         etimef77_(&utime0,&stime0);
      #endif
    }

   /*** Initialization of the algorithms ***/

   ida1=0;
   ida2=1;
   ixlast=0;
   level0=0;
   if(idELL) ndiff=1;

   #if defined (_IBMR2) || defined (AIX) || defined (hpux)
      setcommons(&ndiff,&nonlin,&nreact,&nbc,&ntrue,&idELL,&idpsrc,
                 at,bt,ax,bx,ay,by,wavet,wavex,wavey);
      for(imode=1;imode<=2;imode++)
         matrix3(&imode,&nx,&ny,&level0,ax,bx,ay,by,A[imode-1],ierr);
      initial(&nx,&ny,&idpsrc,&level0,&id_hsrc,ax,bx,ay,by,at,x[0]);
   #else
      setcommons_(&ndiff,&nonlin,&nreact,&nbc,&ntrue,&idELL,&idpsrc,
                 at,bt,ax,bx,ay,by,wavet,wavex,wavey);
      for(imode=1;imode<=2;imode++)
         matrix3_(&imode,&nx,&ny,&level0,ax,bx,ay,by,A[imode-1],ierr);
      initial_(&nx,&ny,&idpsrc,&level0,&id_hsrc,ax,bx,ay,by,at,x[0]);
   #endif

   if(id_hsrc)
      getSRC(nx,ny,level0,interpol,n1,n2,
             ax,bx,ay,by,at,tsbgn,tsend,HEATSRC,src[0],ierr);
   else{
      #if defined (_IBMR2) || defined (AIX) || defined (hpux)
         getrhs(&nx,&ny,&level0,ax,bx,ay,by,at,src[0],ierr);
      #else
         getrhs_(&nx,&ny,&level0,ax,bx,ay,by,at,src[0],ierr);
      #endif
    }


   /*============================================================*/
   if(idELL)
   /*============================================================*/

    { if(idsol!=3 && idsol!=4)
         { printf("heat_trans.c: No case!\n"); *ierr=1; return; }

      printf(" **************************************\n");
      printf(" * Elliptic Problem with Dirichlet BC *\n");
      printf(" * idELL= %d  ndiff=%d                  *\n",idELL,ndiff);
      printf(" **************************************\n");

      sol_ADI_ELL(nt,nx,ny,itmax,idsol,level,maxADI,maxGS,idSGS,
          ndiff,nreact,nbc,idpsrc,n1,n2,interpol,id_hsrc,&ixlast,
          at,bt,ax,bx,ay,by,tol,tsbgn,tsend,eta,
          A5,A3[0],A3[nx*ny*3],A[ida1],A[ida2],x,b,src,wksp,ws,HEATSRC,ierr);
    }

   /*============================================================*/
   else{
   /*============================================================*/

   /*** Call different algorithms ***/

   if(idsol==1)

      sol_LU(nt,nx,ny,itmax,idsol,id_nl,level,
          ndiff,nreact,nbc,idpsrc,n1,n2,interpol,id_hsrc,
          at,bt,ax,bx,ay,by,tol,tsbgn,tsend,
          K4,A5,A[ida1],A[ida2],x,b,src,ws,HEATSRC,&ixlast,ierr);

   else if(idsol==2)

      sol_PCGILU(nt,nx,ny,itmax,idsol,id_nl,idpcond,level,
          ndiff,nreact,nbc,idpsrc,n1,n2,interpol,id_hsrc,idsymm,
          at,bt,ax,bx,ay,by,tol,tsbgn,tsend,
          K4,A5,A[ida1],A[ida2],x,b,src,wksp,HEATSRC,&ixlast,ierr);

   else if(idsol==3)

      #if 1

      sol_ADI(nt,nx,ny,itmax,idsol,level,
          ndiff,nreact,nbc,idpsrc,n1,n2,interpol,id_hsrc,
          at,bt,ax,bx,ay,by,tol,tsbgn,tsend,
          A3[0],A3[nx*ny*3],A[ida1],A[ida2],x,b,src,wksp,HEATSRC,ierr);

      #else

      sol_ADI2(nt,nx,ny,itmax,idsol,level,
          ndiff,nreact,nbc,idpsrc,n1,n2,interpol,id_hsrc,
          at,bt,ax,bx,ay,by,tol,tsbgn,tsend,
          A3[0],A3[nx*ny*3],A[ida1],A[ida2],x,b,src,wksp,ws,HEATSRC,ierr);

      #endif

   else if(idsol==4)

      sol_ADI_II(nt,nx,ny,itmax,idsol,level,maxADI,maxGS,idSGS,
          ndiff,nreact,nbc,idpsrc,n1,n2,interpol,id_hsrc,&ixlast,
          at,bt,ax,bx,ay,by,tol,tsbgn,tsend,eta,
          A5,A3[0],A3[nx*ny*3],A[ida1],A[ida2],x,b,src,wksp,ws,HEATSRC,ierr);

   else if(idsol==5)

      sol_ADI_II_SOR(nt,nx,ny,itmax,idsol,level,maxADI,maxGS,idSGS,
          ndiff,nreact,nbc,idpsrc,n1,n2,interpol,id_hsrc,&ixlast,
          at,bt,ax,bx,ay,by,tol,tsbgn,tsend,eta,
          A5,A3[0],A3[nx*ny*3],A[ida1],A[ida2],x,b,src,wksp,ws,HEATSRC,ierr);

   else{

      if(level>=1) printf(" No solver is given for idsol=%d\n",idsol);
      *ierr=1; return;

    }

   /*============================================================*/
    }
   /*============================================================*/


   /*** Check the computation time ***/

   if(level>=1){
      #if defined (_IBMR2) || defined (AIX) || defined (hpux)
         etimef77(&utime1,&stime1);
      #else
         etimef77_(&utime1,&stime1);
      #endif
      printf(" ElapsedTime=%.3gu %.3gs\n",utime1-utime0,stime1-stime0);
    }


   /*** Numerical analysis: check error @ t=bt ***/

   if(idpsrc==0 && level>=1){

      if(idELL) *bt=0.0;

      #if 0
        numeranal(nt,nx,ny,level,at,bt,ax,bx,ay,by,x[ixlast],ierr);
      #else
        #if defined (_IBMR2) || defined (AIX) || defined (hpux)
           numeranalf(&nt,&nx,&ny,&level,at,bt,ax,bx,ay,by,x[ixlast],ierr);
        #else
           numeranalf_(&nt,&nx,&ny,&level,at,bt,ax,bx,ay,by,x[ixlast],ierr);
        #endif
      #endif
   }

}