sol_adi_ell.c

二维热能方程的求解

#include <stdio.h>
#include <math.h>
#define max(a,b) (a>b)?a:b

typedef  double ATYPE;
typedef  double VTYPE;
typedef  float DTYPE;


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

int   nt,nx,ny,itmax,idsol,level,maxADI,maxGS,idSGS;
int   ndiff,nreact,nbc,idpsrc,n1,n2,interpol,id_hsrc,*ixlast;
VTYPE *at,*bt,*ax,*bx,*ay,*by,*tol,*tsbgn,*tsend,*eta;
ATYPE A5[][nx][5];
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,rho,tworho,pi,ram1,ram2,rhat,optrho;
   VTYPE coefI,coefA1,coefA2,zero,half,two;
   VTYPE top8,bot8,norm8;

   if(level>=1)
      printf("SOL_ADI_ELL: idsol=%d eta=%g\n",idsol,*eta);
   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;

   pi=acos(-1.0);
   ram1=2.*(1.-cos(pi/(nx-1)));
   ram2=2.*(1.-sin(pi/(nx-1)));
   rhat=sqrt(ram1*ram2);
   optrho=hx*hx/rhat;

   rho=optrho;  /****** setting rho *****/

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

   coefI=1.0;
   coefA1=rho;
   coefA2=rho;
   tworho=2.*rho;

  /**********************/
  /* Zero initial guess */
  /**********************/

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


  /**********************/
  /* Scale src by 2.rho */
  /**********************/

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


  /********************/
  /*** Make H and V ***/
  /********************/

   for(iy=0;iy<ny;iy++)
   for(ix=0;ix<nx;ix++){
      H[iy][ix][0]=  -rho*A1[iy][ix][0];
      H[iy][ix][1]=1.-rho*A1[iy][ix][1];
      H[iy][ix][2]=  -rho*A1[iy][ix][2];
     }
   for(ix=0;ix<nx;ix++)
   for(iy=0;iy<ny;iy++){
      V[ix][iy][0]=  -rho*A2[ix][iy][0];
      V[ix][iy][1]=1.-rho*A2[ix][iy][1];
      V[ix][iy][2]=  -rho*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_ADI_ELL: ADI-marching\n");
   n=0; norm8=10.;

   while(norm8>(*tol)){

      n++;

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

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

         if(n==1 || idsol==3){

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

         }
         else{
            for(iy=0;iy<ny;iy++)
            for(ix=0;ix<nx;ix++)
                ws[0][iy][ix]=x[nd0][iy][ix]
                             +(*eta)*(x[nd0][iy][ix]-x[nd1][iy][ix]);

                /*
                ws[0][iy][ix]=two*x[nd0][iy][ix]-x[nd1][iy][ix];
                */

            eval_Ax(nx,ny,level0,zero,zero,coefA2,A1,A2,ws[0],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[1][iy][ix]=src[0][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[1][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

         /*** 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];

         /*** check convergence ***/

         top8=bot8=0.;

         for(iy=1;iy<ny;iy+=2)
         for(ix=1;ix<nx;ix+=2){
            top8=max(top8,fabs(x[nd1][iy][ix]-x[nd0][iy][ix]));
            bot8=max(bot8,fabs(x[nd1][iy][ix]));
          }
         norm8=top8/bot8;

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

    *ixlast=nd1;

    if(level>=1){
        printf(" nx=%d ny=%d rho=%.3g  optrho=%.3g rhat=%.3g\n",
               nx,ny,rho,optrho,rhat);
        printf(" Iter=%d Rel_L8_Norm=%.2g tol=%.2g\n",n,norm8,*tol);
    }


}