sx07.c

解决很多地下水数值模拟问题

 /*PRACTICE 07 (DF-IM-D2SOR)*/
#include <stdio.h>
#include <math.h>
#include<conio.h>
#include<stdlib.h>
void main()
{
	FILE *fp1,*fp2,*fp3,*fp4,*fp5,*fp6,*fp7,*fp8,*fp9,*fp10,*fp11;
	float dx,dy,dt,dtdat,omega,time,qbqli,mxy1t,df,hgs,
		  dfmax,oldval,qwli,hb1li,d;
    float nn[25][12],kmd[3],km[14][6],qbq[18][2],hba[14][2],
		  h[14][6],qw[10][2],qvh[25][12],ww[25][12],mued[2],ee[25][12],
		  ss[25][12],cc[25][12],ff[25][12],mu[14][6];
	int   mdmue[14][6],ijw[10][2],ijba[14][2],mdkm[14][6],ijbq[18][2];
	int   ni,nj,ndkm,ndmue,npba,npbq,nw,ndat,ndt,ni01,nj01,i,j,
	      l,icount,iprint,mdt,mdkmij,mm,mi,mj,mdmuij;
	
        ni=13;
        nj=5;
        ndkm=2;
        ndmue=1;
        npba=13;
        npbq=17;
        nw=1;
        ndat=1;
        ndt=11;
        ni01=ni-1;
        nj01=nj-1;
        dx=100.;
        dy=100.;
        dt=0.5;
        dtdat=100.;
        omega=1.4;
		if((fp1=fopen("Mdmue.sx7","r"))==NULL)
             printf("cant open mdmue.sx7");
		for(i=0;i<=ni;i++)
		{
			for(j=0;j<=nj;j++)
			{   
				fscanf(fp1,"%d",&mdmue[i][j]);
			}
		}
        fclose(fp1);
        fp2=fopen("Mued.sx7","r+");
     	for(i=0;i<=ndmue;i++)
		{
			fscanf(fp2,"%f",&mued[i]);
		}
		fclose(fp2);
		fp3=fopen("Mdkm.sx7","r+");
		for(j=0;j<=nj;j++)
		{
			for(i=0;i<ni;i++)
			{
				fscanf(fp3,"%d",&mdkm[i][j]);
			}
		}
		fclose(fp3);
		fp4=fopen("Kmd.sx7","r+");
		for(i=0;i<=ndkm;i++)
		{     
			fscanf(fp4,"%f",&kmd[i]);
		}
		fclose(fp4);
		fp5=fopen("ijba.sx7","r+");
		for(l=0;l<=npba;l++)
		{
			for(j=0;j<=1;j++)
			{
				fscanf(fp5,"%d",&ijba[l][j]);
			}
		}
		fclose(fp5);
        fp6=fopen("hba.sx7","r");
        for(l=0;l<=npba;l++)
		{
			for(j=0;j<=ndt;j++)
			{
				fscanf(fp6,"%f",&hba[l][j]);
			}
		}
		fclose(fp6);
		/*READ(1,*)((HB1(L,J),J=1,NDAT),L=1,NPB1)
        CLOSE(1)*/
        fp7=fopen("h.sx7","r");
        for(i=0;i<=ni;i++)
		{
			for(j=0;j<=nj;j++)
			{
				fscanf(fp7,"%f",&h[i][j]);
			}
		}
		fclose(fp7);
        /*READ(1,*)((H(I,J),J=1,NJ),I=1,NI)
        CLOSE(1)*/
		fp8=fopen("ijbq.sx7","r");
		for(l=0;l<=npbq;l++)
		{
			for(j=0;j<=1;j++)
			{
				fscanf(fp8,"%d",&ijbq[l][j]);
			}
		}/**/
		fclose(fp8);
        fp9=fopen("qbq.sx7","r");
		for(l=0;l<=npbq;l++)
		{
			for(j=0;j<=ndat;j++)
			{
				fscanf(fp9,"%f",&qbq[l][j]);
			}
		}
		fclose(fp9);
        /*READ(1,*)((QBQ(L,J),J=1,NDAT),L=1,NPBQ)
        CLOSE(1)*/
        fp10=fopen("ijw.sx7","r");
		for(l=0;l<=nw;l++)
		{
			for(j=0;j<=1;j++)
			{
				fscanf(fp10,"%d",&ijw[l][j]);
			}
		}
		fclose(fp10);
        /*READ(1,*)((IJW(L,J),J=1,2),L=1,NW)
        CLOSE(1)*/
        fp11=fopen("qw.SX7","r");
		for(l=0;l<=nw;l++)
		{
			for(j=0;j<=ndat;j++)
			{
				fscanf(fp11,"%f",&qw[l][j]);
			}
		}
		fclose(fp11);

		for(i=0;i<=ni;i++)/*DO 40 I=1,NI*/
        { 
			for(j=0;j<=nj;j++)/*DO 40 J=1,NJ*/
			{
				mdmuij=mdmue[i][j];
                mu[i][j]=mued[mdmuij];
			}
		}
		for(i=0;i<=ni;i++)/*DO 50 I=1,NI*/
		{
			for(j=0;j<=nj;j++)
			{
             mdkmij=mdkm[i][j];
			 km[i][j]=kmd[mdkmij];
			}
		}
    	for(j=1;j<=nj01;j++)/*O 60 J=2,NJ01*/
		{
			for(i=1;i<=ni01;i++)
			{
        ww[i][j]=2*km[i-1][j]*km[i][j]/(km[i-1][j]+km[i][j])*(dy/dx);
        ee[i][j]=2*km[i+1][j]*km[i][j]/(km[i+1][j]+km[i][j])*(dy/dx);
        nn[i][j]=2*km[i][j-1]*km[i][j]/(km[i][j-1]+km[i][j])*(dx/dy);
        ss[i][j]=2*km[i][j+1]*km[i][j]/(km[i][j+1]+km[i][j])*(dx/dy);
			}
		}
/*C       BLOCK 3*/

        icount=1;
        iprint=4;
        time=0.;
        for(mdt=1;mdt<=ndt;mdt++)/*DO 160 MDT=1,NDT*/
		{
			time=time+dt;

        for(i=0;i<=ni;i++)/*DO 65 I=1,NI*/
        {
			for(j=0;j<=nj;j++)/*DO 65 J=1,NJ*/
			{
				qvh[i][j]=0.;
			}
		}/*65      QVH(I,J)=0.0*/
        mm=(int)((time/dtdat)+1);
        d=time-(mm-1)*dtdat;
        for(l=1;l<=npba;l++)/*O 70 L=1,NPB1*/
        {
			hb1li=(hba[l][mm+1]-hba[l][mm])*d/dtdat
             +hba[l][mm];
        mi=ijba[l][0];
        mj=ijba[l][1];
        h[mi][mj]=hb1li;
		}/*70      CONTINUE*/
        for(l=0;l<=npbq;l++)/*DO 90 L=1,NPBQ*/
		{
			qbqli=(qbq[l][mm+1]-qbq[l][mm])*d/dtdat+qbq[l][mm];
			mi=ijbq[l][0];
            mj=ijbq[l][1];
         qvh[mi][mj]=qvh[mi][mj]+qbqli;
		}/*90      CONTINUE*/
       	for(l=0;l<=nw;l++)/*DO 100 L=1,NW*/
		{
			qwli=(qw[l][mm+1]-qw[l][mm])*d/dtdat+qw[l][mm];
			mi=ijw[l][0];
			mj=ijw[l][1];
            qvh[mi][mj]=qvh[mi][mj]+qwli;
		}/*100     CONTINUE*/
        for(i=1;i<=ni01;i++)/*DO 110 I=2,NI01*/
		{
			for(j=1;j<=nj01;j++)
			{
				mxy1t=mu[i][j]*dx*dy/dt; /*DO 110 J=2,NJ01*/
				cc[i][j]=-mxy1t-ww[i][j]-ee[i][j]-nn[i][j]-ss[i][j];
        ff[i][j]=-mxy1t*h[i][j]-qvh[i][j];
			}
		}/*110     CONTINUE*/

/*C       BLOCK 4*/

do
{      dfmax=0.0;
        for(i=1;i<=ni01;i++)
		{
			for(j=1;j<=nj01;j++)
			{
        oldval=h[i][j];
        hgs=(ff[i][j]-ww[i][j]*h[i-1][j]-ee[i][j]*h[i+1][j]-nn[i][j]
			*h[i][j-1]-ss[i][j]*h[i][j+1])/cc[i][j];
        h[i][j]=(1-omega)*oldval+omega*hgs;
        df=fabs(h[i][j]-oldval);
        if(df>dfmax)dfmax=df;
		    }
		}
}while(dfmax>0.001) ;
        if(icount==iprint)
		{
        printf("TIME=%f\n",time);
        for(j=0;j<=nj01;j++)
		{
			for(i=1;i<=ni01;i++)
			{
             printf("H[I][J]=%f-6.2",h[i][j]);
			}
		}
		}
     	else 
           icount=icount+1;
		}/*160     CONTINUE*/
        return ;
	}