sol_pcgilu.c

二维热能方程的求解

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

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

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

{
   int   ix,iy,k,n,iterPCG,total_iter;
   int   nxny,level0,ida,idailu,nd0,nd1;
   VTYPE tn,ht,hx,hy;
   VTYPE coefI,coefA1,coefA2;
   VTYPE half,halfht;
   VTYPE tn0,t,u,scale,sixth;

   if(level>=1) printf("SOL_PCGILU: idsol=%d idpcond=%d idsymm=%d\n",
               idsol,idpcond,idsymm);
   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);
   nxny=nx*ny;
   halfht=0.5*ht;
   sixth=1./6.;

   /*** make A5[][ny][nx][5] ***/

   ida=0;
   idailu=1;

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

   /*** Symmetrizing the Matrix **/
   half=0.5e0;
   if(idsymm>0){
      for(ix=0;ix<nx;ix++){
         for(k=0;k<5;k++) A5[ida][0][ix][k]*=half;
         for(k=0;k<5;k++) A5[ida][ny-1][ix][k]*=half;
       }
      for(iy=0;iy<ny;iy++){
         for(k=0;k<5;k++) A5[ida][iy][0][k]*=half;
         for(k=0;k<5;k++) A5[ida][iy][nx-1][k]*=half;
       }
    }

   ilu0(nx,ny,level0,A5[ida],A5[idailu],ierr);

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

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

   coefI=1.0;
   coefA1=halfht;
   coefA2=halfht;
   total_iter=0;

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

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

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

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

      /*** add the source term to the result ***/
      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++)
         b[iy][ix]+=halfht*(src[nd0][iy][ix]+src[nd1][iy][ix]);


    /*** add nonlinear term ***/

       if(id_nl){

         if(id_nl==1 || n==1){

            #if defined (_IBMR2) || defined (AIX) || defined (hpux)
               rk4step(&nt,&nx,&ny,&idsol,&level0,&ht,&tn0,&hx,&hy,
                    K4[0],x[nd0],ierr);
            #else
               rk4step_(&nt,&nx,&ny,&idsol,&level0,&ht,&tn0,&hx,&hy,
                    K4[0],x[nd0],ierr);
            #endif

            for(iy=0;iy<ny;iy++)
            for(ix=0;ix<nx;ix++){
               b[iy][ix]+=sixth*(K4[0][iy][ix]
                   +2.*(K4[1][iy][ix]+K4[2][iy][ix])+K4[3][iy][ix]);
             }
          }

         else if(id_nl==2){

            scale=1.0;
            #if defined (_IBMR2) || defined (AIX) || defined (hpux)
               extrapol(&nt,&nx,&ny,&idsol,&level0,&ht,&tn0,&hx,&hy,
                   &scale,b,x[nd0],x[nd1],ierr);
            #else
               extrapol_(&nt,&nx,&ny,&idsol,&level0,&ht,&tn0,&hx,&hy,
                   &scale,b,x[nd0],x[nd1],ierr);
            #endif

          }

       }


    /*** initial guess ***/
      if(n>=2)
         for(iy=0;iy<ny;iy++)
         for(ix=0;ix<nx;ix++)
            x[nd1][iy][ix]=2.e0*x[nd0][iy][ix]-x[nd1][iy][ix];
      else
         for(iy=0;iy<ny;iy++)
         for(ix=0;ix<nx;ix++)
            x[nd1][iy][ix]=x[nd0][iy][ix];


    /*** Symmetrize the vector ***/

      if(idsymm>0){
         for(ix=0;ix<nx;ix++){ b[0][ix]*=half; b[ny-1][ix]*=half; }
         for(iy=0;iy<ny;iy++){ b[iy][0]*=half; b[iy][nx-1]*=half; }
       }

      pcg_ilu0(nx,ny,itmax,idpcond,level0,idpsrc,&iterPCG,
               A5[ida],A5[idailu],x[nd1],b,wksp,tol,ierr);

      total_iter+=iterPCG;
    }

      *ixlast=nd1;

      if(level>=1){
         printf(" Total_PCG_Iteration=%d\n",total_iter);
         printf(" n=(%d,%d,%d):",nt,nx,ny);
         printf(" x[4-Corner&Ctr]=%.3g %.3g %.3g %.3g %.3g\n",
             x[nd1][0][0],x[nd1][0][nx-1],x[nd1][ny-1][0],
             x[nd1][ny-1][nx-1],x[nd1][ny/2][nx/2]);
      }

}