sol_adi2.c

二维热能方程的求解

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

void 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,
            H,V,A1,A2,x,b,src,wksp,ws,HEATSRC,ierr)

int   nt,nx,ny,itmax,idsol,level;
int   ndiff,nreact,nbc,idpsrc,n1,n2,interpol,id_hsrc;
VTYPE *at,*bt,*ax,*bx,*ay,*by,*tol,*tsbgn,*tsend;
ATYPE H[][nx][3];
ATYPE V[][ny][3];
ATYPE A1[][nx][3];
ATYPE A2[][ny][3];
ATYPE x[][ny][nx];
ATYPE b[];
ATYPE src[][ny][nx];
ATYPE wksp[][ny];
ATYPE ws[][ny][nx];
DTYPE HEATSRC[];
int   *ierr;

{
   int   ix,iy,k,n,ione;
   int   nx2,nxny,level0,nd0,nd1,iter;
   VTYPE tn,ht,hx,hy;
   VTYPE coefI,coefA1,coefA2,zero,half,two;
   VTYPE halfht,norm,norm1,tmp;

   if(level>=1)
      printf("sol_ADI2: idsol=%d\n",idsol);
   if(level>=2){
      printf(" nt=%d nx=%d ny=%d\n",nt,nx,ny);
      printf(" at=%g bt=%g ax=%g bx=%g ay=%g by=%g\n",*at,*bt,*ax,*bx,*ay,*by);
      printf(" ndiff=%d nreact=%d nbc=%d idpsrc=%d\n",ndiff,nreact,nbc,idpsrc);
      printf(" itmax=%d tol=%g level=%d\n",itmax,*tol,level);
      printf(" n1=%d n2=%d",n1,n2);
      printf(" interpol=%d id_hsrc=%d\n",interpol,id_hsrc);                  
      printf(" tsbgn=%g tsend=%g\n",*tsbgn,*tsend);
      }

   if(level<3)
      level0=0;
   else
      level0=1;

   ht=(*bt-*at)/nt;
   hx=(*bx-*ax)/(nx-1);
   hy=(*by-*ay)/(ny-1);
   nx2=nx*2;
   nxny=nx*ny;
   halfht=0.5*ht;

   ione=1;
   zero=0.e0;
   half=0.5e0;
   two=2.e0;

   coefI=1.0;
   coefA1=halfht;
   coefA2=halfht;


  /********************/
  /*** make H and V ***/
  /********************/

   for(iy=0;iy<ny;iy++)
   for(ix=0;ix<nx;ix++){
      H[iy][ix][0]=  -halfht*A1[iy][ix][0];
      H[iy][ix][1]=1.-halfht*A1[iy][ix][1];
      H[iy][ix][2]=  -halfht*A1[iy][ix][2];
     }
   for(ix=0;ix<nx;ix++)
   for(iy=0;iy<ny;iy++){
      V[ix][iy][0]=  -halfht*A2[ix][iy][0];
      V[ix][iy][1]=1.-halfht*A2[ix][iy][1];
      V[ix][iy][2]=  -halfht*A2[ix][iy][2];
     }

   /*** LU factorization ***/

   #if defined (_IBMR2) || defined (AIX) || defined (hpux)
      for(iy=0;iy<ny;iy++)
         lufac(&ione,&nx,&level0,H[iy]);
      for(ix=0;ix<nx;ix++)
         lufac(&ione,&ny,&level0,V[ix]);
   #else
      for(iy=0;iy<ny;iy++)
         lufac_(&ione,&nx,&level0,H[iy]);
      for(ix=0;ix<nx;ix++)
         lufac_(&ione,&ny,&level0,V[ix]);
   #endif


   /************************************/
   /*** Now, ready for time-marching ***/
   /************************************/

   if(level>=1) printf("sol_ADI2: time-marching\n");

   for(n=1;n<=nt;n++){

      tn=ht*n;
      if(level>=3) printf(" n=%d tn=%g\n",n,tn);

      nd0=(n+1)%2;
      nd1=n%2;

      /***************/
      /*** get RHS ***/
      /***************/

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

         for(iy=0;iy<ny;iy++)
         for(ix=0;ix<nx;ix++)
            src[nd0][iy][ix]=halfht*(src[nd0][iy][ix]+src[nd1][iy][ix]);

         eval_Ax(nx,ny,level0,zero,zero,coefA2,A1,A2,x[nd0],ws[1],ierr);
         eval_Ax(nx,ny,level0,zero,coefA1,zero,A1,A2,ws[1],ws[0],ierr);

         for(iy=0;iy<ny;iy++)
         for(ix=0;ix<nx;ix++)
            src[nd0][iy][ix]+=ws[0][iy][ix];

      /***********/
      /* X-SWEEP */
      /***********/

         /*** evaluate "coefI*I+coefA1*A1+coefA2*A2" ***/
         eval_Ax(nx,ny,level0,coefI,coefA1,coefA2,A1,A2,x[nd0],x[nd1],ierr);

         for(iy=0;iy<ny;iy++)
         for(ix=0;ix<nx;ix++)
            x[nd1][iy][ix]+=src[nd0][iy][ix];

         /*** forward and back substitutions ***/
         #if defined (_IBMR2) || defined (AIX) || defined (hpux)
            for(iy=0;iy<ny;iy++)
               substit(&ione,&nx,&level0,H[iy],x[nd1][iy]);
         #else
            for(iy=0;iy<ny;iy++)
               substit_(&ione,&nx,&level0,H[iy],x[nd1][iy]);
         #endif

      /***********/
      /* Y-SWEEP */
      /***********/

         /*** Note: the indices are changed for the transpose. ***/

         for(ix=0;ix<nx;ix++)
         for(iy=0;iy<ny;iy++)
            wksp[ix][iy]=x[nd1][iy][ix];

         /*** forward and back substitutions ***/
         #if defined (_IBMR2) || defined (AIX) || defined (hpux)
            for(ix=0;ix<nx;ix++)
               substit(&ione,&ny,&level0,V[ix],wksp[ix]);
         #else
            for(ix=0;ix<nx;ix++)
               substit_(&ione,&ny,&level0,V[ix],wksp[ix]);
         #endif

         /*** Compute the difference: ws[0]=(u^{n+2/3}-u^n).
              Save the solution of the standard ADI.
              Return back to the row-wise ordering ***/

         for(iy=0;iy<ny;iy++)
         for(ix=0;ix<nx;ix++){
            x[nd1][iy][ix]=wksp[ix][iy];
            }

   }
   /****************************/
   /* END OF THE TIME STEPPING */
   /****************************/

   if(level>=1) printf(" n=(%d,%d,%d)\n",nt,nx,ny);

}