一维对流扩散问题TDMA算法

/*
一维对流扩散问题TDMA算法，支持中心差分，迎风，混合，指数，乘方格式，还可
以输出精确解

用法
直接编译运行，不带参数运行时会自动输出使用方法:

upwind –left=(左端点处的值) –right=(右端点处的值) –pdelta=(网格Pe数)
–scheme=(格式，可选center/upwind/mixed/exponent/power/precise) –output=(输出文
件名)

如：
upwind –left=100 –right=200 –pdelta=5 –scheme=upwind –output=Pe5Upwind.dat
参数顺序无关紧要。

注意现在的代码输出相对值(phi-phi0)/(phi_L-phi0)，想输出绝对值phi的自己改代码吧:)
另存为upwind.c编译即可
//*/
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>

double max(double number1,double number2);
void parseParameter(int *argc, char **argv, double *PHI_0, double *PHI_L, double *P_Delta, char *scheme, char *outputfilename);

int main(int argc, char **argv)	
{
	
  const int GRIDAMOUNT=20;//网格数目为20，由于两端点已知，实际未知数个数为19个。
  double PHI_0=100;
	
  double PHI_L=200;
  double P_Delta=1; //网格Pe数，取为1,5,10等。
  char scheme[10]="upwind";
	
  int PRECISE=0;
  char outputfilename[100]="output.dat";
	
  double a_E[GRIDAMOUNT-1];
  double a_W[GRIDAMOUNT-1];
	
  double a_P[GRIDAMOUNT-1];
  double A_P_Delta[GRIDAMOUNT-1];
	
  double P_TDMA[GRIDAMOUNT-1];
  double Q_TDMA[GRIDAMOUNT-1];
	
  double PHI[GRIDAMOUNT-1];
	
  int i;
	
  FILE *dataoutput;
	
  parseParameter(&argc, argv, &PHI_0, &PHI_L, &P_Delta, scheme, outputfilename);
	
  if(!strcmp(scheme,"upwind"))
    {
      for (i=0;i<GRIDAMOUNT-1;i++){A_P_Delta[i]=1;}
	
    }
  else if(!strcmp(scheme,"center"))
	
    {
      for (i=0;i<GRIDAMOUNT-1;i++){A_P_Delta[i]=1-0.5*fabs(P_Delta);}
	
    }
  else if(!strcmp(scheme,"mixed"))
	
    {
      for (i=0;i<GRIDAMOUNT-1;i++){A_P_Delta[i]=max(0,1-0.5*fabs(P_Delta));}
	
    }
  else if(!strcmp(scheme,"exponent"))
	
    {
      for (i=0;i<GRIDAMOUNT-1;i++){A_P_Delta[i]=fabs(P_Delta)/(exp(fabs(P_Delta))-1);}
	
    }
  else if(!strcmp(scheme,"power"))
	
    {
      for (i=0;i<GRIDAMOUNT-1;i++){A_P_Delta[i]=max(0,pow((1-0.1*fabs(P_Delta)),5));}
	
    }
  else if(!strcmp(scheme,"precise"))
	
    {
      PRECISE=1;
    }
  else
    {
	
      fprintf(stderr,"\nUnknown scheme %s!\nAccepted schemes are: upwind|center|mixed|exponent|power|precise\n",scheme);
      exit(2);
	
    }
  //      A_P_Delta[i]=1-0.5*fabs(P_Delta);
	
  //      A_P_Delta[i]=1;
  //      A_P_Delta[i]=max(0,1-0.5*fabs(P_Delta));
	
  //      A_P_Delta[i]=fabs(P_Delta)/(exp(fabs(P_Delta))-1);
	
  //      A_P_Delta[i]=max(0,pow((1-0.1*fabs(P_Delta)),5));
	
  /*
    这里A(|P_Delta|)的取法可见陶书第一版P201，各种格式的取法如下：
    中心差分：1-0.5*fabs(P_{\Delta})
    迎风差分：1
    混合格式：max(0,1-0.5*fabs(P_{\Delta}))
    指数格式：fabs(P_{\Delta})/(exp(fabs(P_{\Delta}))-1)
    乘方格式：max(0,(1-0.1*fabs(P_{\Delta}))^5)
    注意网格P_Delta和整体的P_e的不同！
	
    参考精确解：(exp(P_Delta*(i+1))-1)/(exp(P_Delta*GRIDAMOUNT)-1)*(PHI_L-PHI_0)+PHI_0
  */
	
  /* 下面开始对诸格式均一致的计算过程 */
	
  /* a_W,a_E,a_P赋值 */
  for (i=0;i<GRIDAMOUNT-1;i++)
	
    {
      a_E[i]=A_P_Delta[i]+max(0,-P_Delta);//此处P_Delta在东西界面上相同，所以不区分。
      a_W[i]=A_P_Delta[i]+max(P_Delta,0);
	
      a_P[i]=a_E[i]+a_W[i];//后边的+Fe-Fw项由常物性假定而略去
    }
	
  /* TDMA算法求解这一方程组 */
  P_TDMA[0]=a_E[0]/a_P[0];
	
  Q_TDMA[0]=a_W[0] * PHI_0 / a_P[0];
	
  for (i=1;i<GRIDAMOUNT-2;i++)
	
    {
      P_TDMA[i]=a_E[i]/(a_P[i]-a_W[i]*P_TDMA[i-1]);
	
      Q_TDMA[i]=(a_W[i]*Q_TDMA[i-1])/(a_P[i]-a_W[i]*P_TDMA[i-1]);
	
    }
  P_TDMA[GRIDAMOUNT-2]=a_E[GRIDAMOUNT-2]/(a_P[GRIDAMOUNT-2]-a_W[GRIDAMOUNT-2]*P_TDMA[GRIDAMOUNT-3]);
	
  Q_TDMA[GRIDAMOUNT-2]=(PHI_L*a_E[GRIDAMOUNT-2]+a_W[GRIDAMOUNT-2]*Q_TDMA[GRIDAMOUNT-3])/(a_P[GRIDAMOUNT-2]-a_W[GRIDAMOUNT-2]*P_TDMA[GRIDAMOUNT-3]);
	
  PHI[GRIDAMOUNT-2]=Q_TDMA[GRIDAMOUNT-2];
	
  for (i=GRIDAMOUNT-3;i>=0;i--)
	
    {
      PHI[i]=P_TDMA[i]*PHI[i+1]+Q_TDMA[i];
	
    }
	
  /* 注意：以下为精确解的计算，仅当需要精确解时方才去掉注释语句  */
  if (PRECISE){
    for (i=0;i<GRIDAMOUNT-1;i++)
	
      {
	PHI[i]=(exp(P_Delta*(i+1))-1)/(exp(P_Delta*GRIDAMOUNT)-1)*(PHI_L-PHI_0)+PHI_0;
	
      }
  }
	
  //打印输出，并与精确解比较
  //  printf("\n PHI\tPHI_PRICISE\tP_TDMA\tQ_TDMA"); //标题行
  if((dataoutput=fopen(outputfilename,"wt"))==NULL)
	
    {
      printf("\nerror on open output file!");
      exit(1);
	
    }
  fprintf(dataoutput,"\n %d\t%f",1, (PHI_0-PHI_0)/(PHI_L-PHI_0));//输出相对值 (PHI-PHI_0)/(PHI_L-PHI_0)
	
  for (i=0;i<GRIDAMOUNT-1;i++)
	
    {
      fprintf(dataoutput,"\n %d\t%f", i+2, (PHI[i]-PHI_0)/(PHI_L-PHI_0));
	
    }
  fprintf(dataoutput,"\n %d\t%f\n",GRIDAMOUNT+1,(PHI_L-PHI_0)/(PHI_L-PHI_0));
	
  fclose(dataoutput);
  /*测试strcpy
    printf(outputfilename);
    printf("\n");
    printf(scheme);
    printf("\n");
    strcpy(outputfilename,scheme);
    printf(outputfilename);
    printf("\n");
    printf(scheme);
    printf("\n");
  */
	
  return 0;
}
	
double max (double number1,double number2)
	
{
  double maxvalue=0;
  maxvalue=number2;
  if (number1>=number2)
	
    {
      maxvalue=number1;
    }
  return maxvalue;
	
}
	
void parseParameter(int *argc, char **argv, double *PHI_0, double *PHI_L, double *P_Delta, char *scheme, char *outputfilename)
	
{
  int i=0;
  int pointerofequ=0;
	
  if (*argc<=1)
    {
      fprintf(stderr,"\n Usage:\n%s –scheme SCHEME\n\t –pdelta P_Delta\n\t –left VALUE-OF-LEFT\n\t –right VALUE-OF_RIGHT\n\t –output OUTPUTFILE\n",argv[0]);
	
      exit(0);
    }
  for (i=1;i<*argc;i++)
	
    {
      if ((argv[i][0]=='-')&&(argv[i][1]=='-'))
	
	{
	  argv[i]+=2;
	  pointerofequ=strcspn(argv[i],"=");
	
	  if(!strncmp(argv[i],"pdelta",pointerofequ))
	
	    {
	      argv[i]+=pointerofequ+1;
	      *P_Delta=atol(argv[i]);
	
	    }
	  else if(!strncmp(argv[i],"left",pointerofequ))
	
	    {
	      argv[i]+=pointerofequ+1;
	      *PHI_0=atol(argv[i]);
	
	    }
	  else if(!strncmp(argv[i],"right",pointerofequ))
	
	    {
	      argv[i]+=pointerofequ+1;
	      *PHI_L=atol(argv[i]);
	
	    }
	  else if(!strncmp(argv[i],"scheme",pointerofequ))
	
	    {
	      argv[i]+=pointerofequ+1;
	      strcpy(scheme,argv[i]);
	
	    }
	  else if(!strncmp(argv[i],"output",pointerofequ))
	
	    {
	      argv[i]+=pointerofequ+1;
	      strcpy(outputfilename,argv[i]);
	
	    }
	  else
	    {
	      fprintf(stderr,"\n Unrecognized parameters: –%s\n",argv[i]);
	
	      fprintf(stderr,"\n Usage:\n%s –scheme SCHEME\n\t –pdelta P_Delta\n\t –left VALUE-OF-LEFT\n\t –right VALUE-OF_RIGHT\n\t –output OUTPUTFILE\n",argv[0]);
	
	      exit(0);
	    }
	}
    }
  return;
	
}