sol_lu.c

二维热能方程的求解

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

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

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

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

   if(level>=1) printf("SOL_LU: idsol=%d id_nl=%d\n",idsol,id_nl);
   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.5e0*ht;
   sixth=1./6.;

   /*** make A5[ny][nx][2*nx+1] ***/

   for(iy=0;iy<ny;iy++)
      for(ix=0;ix<nx;ix++)
        for(k=0;k<=nx2;k++)
           A5[iy][ix][k]=0.0;

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

   /*** LU factorization ***/

   #if defined (_IBMR2) || defined (AIX) || defined (hpux)
      lufac(&nx,&nxny,&level0,A5);
   #else
      lufac_(&nx,&nxny,&level0,A5);
   #endif

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

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

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


   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" ***/
      eval_Ax(nx,ny,level0,coefI,coefA1,coefA2,A1,A2,x[nd0],x[nd1],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++)
         x[nd1][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++)
               x[nd1][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.;
            #if defined (_IBMR2) || defined (AIX) || defined (hpux)
               extrapol(&nt,&nx,&ny,&idsol,&level0,&ht,&tn0,&hx,&hy,
                   &scale,x[nd1],x[nd0],ws,ierr);
            #else
               extrapol_(&nt,&nx,&ny,&idsol,&level0,&ht,&tn0,&hx,&hy,
                   &scale,x[nd1],x[nd0],ws,ierr);
            #endif

          }

       }

    /*** save for next step ***/

      if(id_nl==2)
      for(iy=0;iy<ny;iy++)
      for(ix=0;ix<nx;ix++)
         ws[iy][ix]=x[nd0][iy][ix];


    /*** forward and back substitutions ***/
      #if defined (_IBMR2) || defined (AIX) || defined (hpux)
          substit(&nx,&nxny,&level0,A5,x[nd1]);
      #else
          substit_(&nx,&nxny,&level0,A5,x[nd1]);
      #endif


    }

      *ixlast=nd1;

      if(level>=1){
         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]);
      }

}