solve.c

二维SIMPLEC算法

#include "global_var.h" 
void solve()
{
 int jstf,it1,it2,jt1,jt2,n,nt,i,j,istf,ii,jj;
 double bl,blc,blm,blp,denom,temp,deno;     
      printdebuginfo();
      istf=ist-1;
      jstf=jst-1;
      it1=l2+ist;
      it2=l3+ist;
      jt1=m2+jst;
      jt2=m3+jst;
	  
 for (nt=1;nt<=ntimes[nf];nt++)
	 for (n=nf;n<=nf;n++)
	 {
		 if (lblk[nf])
		 {
			 //以列为基本单位的块修正。 
			  pt[istf]=0.;
              qt[istf]=0.;
              for (i=ist;i<=l2;i++)
			  {
			  bl=0.;
              blp=0.;
              blm=0.;
              blc=0.;
              for (j=jst;j<=m2;j++)
			  {
			  bl=bl+ap[i][j];
              if(j!=m2) bl=bl-ajp[i][j];
              if(j!=jst) bl=bl-ajm[i][j];
              blp=blp+aip[i][j];
              blm=blm+aim[i][j];
              blc=blc+con[i][j]+aip[i][j]*f[i+1][j][n]+aim[i][j]*f[i-1][j][n]+ajp[i][j]*f[i][j+1][n]+ajm[i][j]*f[i][j-1][n]-ap[i][j]*f[i][j][n];
           	  }
              denom=bl-pt[i-1]*blm;
              deno=1.e15;
              if(fabs(denom/bl)<1.e-10) 
				  denom=1.e20*deno;
              pt[i]=blp/denom;
              qt[i]=(blc+blm*qt[i-1])/denom;
           
			  } //for i
              bl=0.;
              for (ii=ist;ii<=l2;ii++)
			  {
			  i=it1-ii;
              bl=bl*pt[i]+qt[i];
              for (j=jst;j<=m2;j++)
			     f[i][j][n]=f[i][j][n]+bl;
			  }
              //以行为基本单位的块修正。
              pt[jstf]=0.;
              qt[jstf]=0.;
              for (j=jst;j<=m2;j++)
			  {
              bl=0.;
              blp=0.;
              blm=0.;
              blc=0.;
              for (i=ist;i<=l2;i++)
			  {
              bl=bl+ap[i][j];
              if(i!=l2) bl=bl-aip[i][j];
              if(i!=ist) bl=bl-aim[i][j];
              blp=blp+ajp[i][j];
              blm=blm+ajm[i][j];
              blc=blc+con[i][j]+aip[i][j]*f[i+1][j][n]+aim[i][j]*f[i-1][j][n] +ajp[i][j]*f[i][j+1][n]+ajm[i][j]*f[i][j-1][n]-ap[i][j]*f[i][j][n];
           	  }
              denom=bl-pt[j-1]*blm;
              if (fabs(denom/bl)<1.e-10) denom=1.e20*deno;
              pt[j]=blp/denom;
              qt[j]=(blc+blm*qt[j-1])/denom;
           	  }
              bl=0.;
			  for (jj=jst;jj<=m2;jj++)
			  {
              j=jt1-jj;
              bl=bl*pt[j]+qt[j];
              for (i=ist;i<=l2;i++)
                    f[i][j][n]=f[i][j][n]+bl;
			  }
		 }  //if lblk(nf) 块修正结束
     printdebuginfo();
    
	 //第一次TDMA求解。求解顺序："外循环=从上到下(j)"
	 for (j=jst;j<=m2;j++)
	  {
      pt[istf]=0.;
      qt[istf]=f[istf][j][n];
      for (i=ist;i<=l2;i++)
	  {
      denom=ap[i][j]-pt[i-1]*aim[i][j];
      pt[i]=aip[i][j]/denom;
      temp=con[i][j]+ajp[i][j]*f[i][j+1][n]+ajm[i][j]*f[i][j-1][n];
      qt[i]=(temp+aim[i][j]*qt[i-1])/denom;
	  }
	  for (ii=ist;ii<=l2;ii++)
	  {
       i=it1-ii;
       f[i][j][n]=f[i+1][j][n]*pt[i]+qt[i];
	  }
	}  //for j
     printdebuginfo();
	 
	 //第二次TDMA求解。求解顺序："外循环=从下到上(j)"
     for (jj=jst;jj<=m3;jj++)
	 {
      j=jt2-jj;
      pt[istf]=0.;
      qt[istf]=f[istf][j][n];
      for (i=ist;i<=l2;i++)
	  {
	  denom=ap[i][j]-pt[i-1]*aim[i][j];
      pt[i]=aip[i][j]/denom;
      temp=con[i][j]+ajp[i][j]*f[i][j+1][n]+ajm[i][j]*f[i][j-1][n];
      qt[i]=(temp+aim[i][j]*qt[i-1])/denom;
	  }
      for (ii=ist;ii<=l2;ii++)
	  {
	   i=it1-ii;
       f[i][j][n]=f[i+1][j][n]*pt[i]+qt[i];
	  }
	 } //for jj
     printdebuginfo();
	 
	 //第三次TDMA求解。求解顺序："外循环=从左到右(i)"
	 for (i=ist;i<=l2;i++)
	 {
      pt[jstf]=0.;
      qt[jstf]=f[i][jstf][n];
      for (j=jst;j<=m2;j++)
	  {
	   denom=ap[i][j]-pt[j-1]*ajm[i][j];
      pt[j]=ajp[i][j]/denom;
      temp=con[i][j]+aip[i][j]*f[i+1][j][n]+aim[i][j]*f[i-1][j][n];
      qt[j]=(temp+ajm[i][j]*qt[j-1])/denom;
	  }
      for (jj=jst;jj<=m2;jj++)
	  {
	    j=jt1-jj;
        f[i][j][n]=f[i][j+1][n]*pt[j]+qt[j];
	  }
	 }  //for i
     printdebuginfo();

	 //第四次TDMA求解。求解顺序："外循环=从右到左(i)"
	 for  (ii=ist-1;ii<=l3;ii++)
	 {
      i=it2-ii;
      pt[jstf]=0.;
      qt[jstf]=f[i][jstf][n];
      for (j=jst;j<=m2;j++)
	  {
      denom=ap[i][j]-pt[j-1]*ajm[i][j];
      pt[j]=ajp[i][j]/denom;
      temp=con[i][j]+aip[i][j]*f[i+1][j][n]+aim[i][j]*f[i-1][j][n];
      qt[j]=(temp+ajm[i][j]*qt[j-1])/denom;
	  }
     for (jj=jst;jj<=m2;jj++)
	 {
      j=jt1-jj;
      f[i][j][n]=f[i][j+1][n]*pt[j]+qt[j];
	 }
	 }  
     printdebuginfo(); //第四次TDMA求解。

 } //for nt
	 for (j=2;j<=m2;j++)
		 for (i=2;i<=l2;i++)
		 {
           con[i][j]=0.;  //sc=0
           ap[i][j]=0.;   //sp=0
		 }
		 
}