nl.c

一个关于多重网格的程序

  free(v);
  free(d);
  free(dc);
  free(uc);
  free(vc);

	     
} /* mg */ 


void ic(double *u,int n,int t) // initial condition
{
  int i,j;
  double x,y,h,timebefore;

  h=(xmax-xmin)/(double) n;

  timebefore=(double)(t-1)*tau;
  for(i=0;i<n;i++){
    x=i*h+xmin; 
    for(j=0;j<n;j++){
      y=j*h+xmin; 
      u[index2d(i,j,n)]=sin(x)*sin(y)*sin(timebefore);
    }//i
  }//j
}/* ic */

double error(double *u,double *uback, int n,int *imx,int *jmx)
{   double diff,diffmx,ave;
 int i,j,ic,k;		/* calculates the difference between u now &
			   u before, as one way to calculate error */
 ave=0;
 diffmx=0;
 ic=0;
 for(i=0;i<n;i++){
   for(j=0;j<n;j++){
     diff=u[index2d(i,j,n)]-uback[index2d(i,j,n)];
     if(diff < 0) diff=-diff;
     ave+=diff;
     ic+=1;
     if(diff > diffmx){
       diffmx=diff; 	
       *imx=i;
       *jmx=j;
     } //if
   }//j
 }//i
 ave=ave/ic;
 return (diffmx);
}/*error*/


void pooofle(double *u,int n,double *f,double *nlin)
{
  double *v,*uback,*dc,diff,tol,*nlinc;
  double *d,*vc,*uc,*fc,sumu,sumu2;
  int i,j,ir,imx,jmx,nc;
  int nmem,ncmem,nfmem,count,w,nsq;
  FILE *f1;

  nmem=(n+1)*(n+1);
  nsq=n*n;

  for(i=0;i<nmem;i++){
    u[i]=0.;
  }

  nc=n/2;

  if(n > nmin){

    for(i=0;i<nrep1;i++){
      newton(nlin,n,f,u); // note that it is u that changes
    }       

    d=malloc(nmem*sizeof(double));
    if(d==NULL)printf(" oopsies can't allocate d %d\n",__LINE__);
    for( i=0;i<nmem;i++){ 
      d[i] =0.;
    }

    defectnl(d,u,n,f);
    ncmem=(nc+1)*(nc+1);

    dc=malloc(ncmem*sizeof(double));
    if(dc==NULL)printf(" oopsies can't allocate dc %d\n",__LINE__);
    nlinc=malloc(ncmem*sizeof(double));
    if(nlinc==NULL)printf(" oopsies can't allocate nlinc %d\n",__LINE__);
    vc=malloc(ncmem*sizeof(double));
    if(vc==NULL)printf(" oopsies can't allocate vc %d\n",__LINE__);
    for( i=0;i<ncmem;i++){ 
      dc[i] =0.;
      nlinc[i] =0.;
      vc[i] =0.;
    }


    restrct(dc,d,nc); /* dc is defect on coarse grid, just calculated it */
    restrct2(nlinc,nlin,nc);
    /* restricts the original function u to a coarser grid, since in the
     * linearized version of things need u in the operator K 
     * Here the original function u is called 'nlin'. confusing change
     * of notation
     */

    pooofle(vc,nc,dc,nlinc);

    free(dc);
    free(d);
    free(nlinc);

    v=malloc(nmem*sizeof(double));
    if(v==NULL) printf(" oopsies can't allocate v  %d\n",__LINE__);
    for( i=0;i<nmem;i++){ 
      v[i] =0.;
    }

    interpolate(vc,v,nc,n); /*  interpolate vc up to the fine grid */
    free(vc);

    for(i=0;i<nmem;i++){
      u[i]+=v[i];
    }

    free(v);

    for(i=0;i<nrep2;i++){ 	// postsmoothing
      newton(nlin,n,f,u);
    }
       
  }else{ /* ie do this when at the coarsest grid*/

    //      tol=10.e-10;
    //diff=5*tol;
    //ir=0;

    uback=malloc(nmem*sizeof(double));
    if(uback==NULL)printf(" oopsies can't allocate uback %d\n",__LINE__);
    //         while(diff > tol){    
    /* this is usually here in MG codes
     * it's the 'exact solution on coarse grid' bit
     * I've just found it makes no difference to 
     * the answer in the current code & it takes time
     *
     */
    diff=0;		
    ir=ir++;		
    sumu=0.;
    sumu2=0.;
    for(i=0;i<nsq;i++){		// need to get the avg of u to be 0 
      uback[i]=u[i];		// this is because of the periodic bc's
      sumu +=u[i];		// this acts as a 2nd boundary condition
    }/* i */
    sumu=sumu/nsq;
    for(i=0;i<nsq;i++){
      u[i] -= sumu ;
      //   sumu2 +=u[i];  // sum2 was just a check on u
    }/* i */

    newton(nlin,n,f,u);
    diff=error(u,uback,n,&imx,&jmx);

    //         }/* while */
    free(uback);
  }/* end if/else */
}/* end of poofle */


int index2d(int i,int  j,int  stride)	// this is the wrap-around condition for
{   					// the periodic boundaries
  if(i<0) i+=stride;
  else if(i>=stride) i-=stride;
  if(j<0) j+=stride;
  else if(j>=stride) j-=stride;
  return(i+stride*j);
}

void interpolate(double *uc,double *uf,int nc,int nf)
{
  int i1,i,j,j1,n;
  int ic,jc;
  /*       first equate the matching points on grid  */

  n=nc;
  for( ic=0;ic<nc;ic++){  /* remember uc,nc is coarse grid, uf nf is fine */
    for( jc=0;jc<nc;jc++){       
      i1=ic*2;
      j1=jc*2;
      uf[index2d(i1,j1,nf)]=uc[index2d(ic,jc,nc)]; /* for uf this is even i, even j */
    }//j1
  }//i1

  /* now match up i pts between defined from previous pts - these are 'odd i', even j */
  n=nf;
  for( i=0;i<n;i+=2){ /* i should increment by 2 */
    for( j=1;j<n;j+=2){       
      uf[index2d(i,j,nf)]=(uf[index2d(i,j-1,nf)] + uf[index2d(i,j+1,nf)])*0.5;
    }//j
  }//i

  for( i=1;i<n;i+=2){
    for( j=0;j<n;j++){
      uf[index2d(i,j,nf)]=(uf[index2d(i-1,j,nf)] + uf[index2d(i+1,j,nf)])*0.5; /** note that i,j has been reversed **/
    }//i
  }//j

}/* interpolate */

void restrct(double *uc, double *uf, int n)
{
  int i, j, i2, j2,nf;
  nf=2*n;
      
  /*first equate the matching points on grid */
  /*uc=coarse, uf=fine */
  for( i=0;i<n;i++){ // n is the # of pts on the coarse grid 
    for( j=0;j<n;j++) {
      i2=2*i;
      j2=2*j;
      uc[index2d(i,j,n)]=uf[index2d(i2,j2,nf)] *4.
	+2.*uf[index2d(i2,j2+1,nf)] +2.*uf[index2d(i2+1,j2,nf)]
	+2.*uf[index2d(i2-1,j2,nf)]+2.*uf[index2d(i2,j2-1,nf)]
	+uf[index2d(i2-1,j2+1,nf)] +uf[index2d(i2+1,j2-1,nf)]
	+uf[index2d(i2-1,j2-1,nf)] +uf[index2d(i2+1,j2+1,nf)];
      uc[index2d(i,j,n)]=uc[index2d(i,j,n)]/16.;

      // alternate schemes:
      //           uc[index2d(i,j,n)]=uf[index2d(i2,j2,nf)] ;
      //           or
      //           uc[index2d(i,j,n)]=uf[index2d(i2,j2,nf)] *4.
      //     +uf[index2d(i2,j2+1,nf)] +uf[index2d(i2+1,j2,nf)]
      //     +uf[index2d(i2-1,j2,nf)]+uf[index2d(i2,j2-1,nf)]
      //     uc[index2d(i,j,n)]=uc[index2d(i,j,n)]/8;
		    

    }//j
  }//i

} /* end  restrct */

void restrct2(double *uc, double *uf, int n)
{
  int i, j, i2, j2,nf;
  nf=2*n;
      
  for( i=0;i<n;i++){ /* n is the # of pts on the coarse grid */
    for( j=0;j<n;j++) {
      i2=2*i;
      j2=2*j;
      uc[index2d(i,j,n)]=uf[index2d(i2,j2,nf)]*4.
	+uf[index2d(i2,j2+1,nf)] +uf[index2d(i2+1,j2,nf)]
	+uf[index2d(i2-1,j2,nf)]+uf[index2d(i2,j2-1,nf)];
      uc[index2d(i,j,n)]=uc[index2d(i,j,n)]/8.;

    }/* loop over j */
  }/* loop over i */

} /* end  restrct2 */

double errordefect( double *d, int n) // calculates one norm of the 
	                             // matrix d
{
  int i,j;
  double diff;

  diff=0.;
  for(i=0; i<n; i++){
    for(j=0; j<n; j++){
      diff += d[index2d(i,j,n)]*d[index2d(i,j,n)];
    }
  }
  return(diff);
}/* errordefect */

void newton( double *u, int n, double *f, double *v)
{
  /* linearized newtons method for nonlinear pde  
   * see p. 149 trottenberg for example */ 

  int i,j,nmem;
  double h,h2,d2;
  double nonlin,denom,*v2;

  h=(xmax-xmin)/(double) n; 
  //this is h : h= (xmax-xmin)/n, assuming xmax/min=ymax/min
  h2=h*h;     //  if you change this, you must also change d2 in defect !!!
  d2=1./h2;

  for(i=0;i<n;i++){
    for(j=0;j<n;j++){
      denom= -4.*d2 +invtau +df(u[index2d(i,j,n)]);
      v[index2d(i,j,n)]=
	(f[index2d(i,j,n)]-(v[index2d(i-1,j,n)]+v[index2d(i,j-1,n)]
       + v[index2d(i+1,j,n)]+v[index2d(i,j+1,n)] )*d2 )/denom;
    }
  }

}/*newton*/


void defectnl(double *d, double *u,int n,double *f) 
{
  int i,j;
  double h,h2,d2,temp,temp1,temp2,temp3,b1;

  h=(xmax-xmin)/(double) n; //this is h : h= (xmax-xmin)/n
  h2=h*h;
  d2=1./h2;
  for( i=0;i<n;i++){ 
    for( j=0;j<n;j++){  // d = f - Nu
      d[index2d(i,j,n)]=f[index2d(i,j,n)] -
	+ (  (u[index2d(i-1,j,n)]+u[index2d(i,j-1,n)]
	      +u[index2d(i+1,j,n)]+u[index2d(i,j+1,n)] -4*u[index2d(i,j,n)])*d2
	     +invtau*u[index2d(i,j,n)]
	     +fn(u[index2d(i,j,n)]) );

    }
  }
}/* defectnl */

void funct(double *f,double *u,int n, int t)
{
  int i,j,k;
  double x,y,z,h,t1;
  double d2,sxsy,sint,sint1,uij,fij,fm1,fm,cost,cost1;
  t1=(t-1)*tau;
  d2=(double)n/(xmax-xmin);
  d2=d2*d2;
  h=(xmax-xmin)/(double) n;
  for(i=0;i<n;i++){
    x=i*h +xmin;
    for(j=0;j<n;j++){ 
      y=j*h +xmin;
      sxsy=sin(x)*sin(y); // definitions to make things easier
      sint=sin(t*tau);
      sint1=sin(t1);
      cost=cos(t*tau);
      cost1=cos(t1);
      uij=u[index2d(i,j,n)];
      // define the functions f^m+1, f^m here to makes things easier
      // eqn is : Nu +f = d_t u  => f=d_t u-Nu
      // need f at this time (t), and previous time (t-1)
      // see comments at begining for more explanation
      fm  = sxsy*cost   + 2*sxsy*sint  - fn(sxsy*sint);
      fm1 = sxsy*cost1 + 2*sxsy*sint1 - fn(sxsy*sint1);
      fij=(u[index2d(i+1,j,n)]+u[index2d(i-1,j,n)]+
	   u[index2d(i,j+1,n)]+u[index2d(i,j-1,n)] -4.*uij)*d2;
      fij+=uij*(2./tau) + fn(uij) + fm1 +fm;
      f[index2d(i,j,n)]=-fij;
  
    }//j
  }//i
} /* funct */ 

void newtonnl( double *u, int n, double *f)
{ //newtons method directly for pde
  int i,j,nmem; double h,h2,d2,*u2;
  double nonlin,denom;

  h=(xmax-xmin)/(double) n;
  h2=h*h;     //  if you change this, you must also change d2 in defect !!!
  d2=1./h2;
  for(i=0;i<n;i++){ 
    for(j=0;j<n;j++){ 
      denom= -4.*d2 +invtau+df(u[index2d(i,j,n)]);
      nonlin=u[index2d(i,j,n)]*invtau +fn(u[index2d(i,j,n)]) -f[index2d(i,j,n)];
      u[index2d(i,j,n)]=
	u[index2d(i,j,n)]
	-( (u[index2d(i-1,j,n)]+u[index2d(i,j-1,n)]+ u[index2d(i+1,j,n)]
	   +u[index2d(i,j+1,n)] -4*u[index2d(i,j,n)])*d2 +nonlin )/denom;
    }//j
  }//i

}/*newtonnl*/