nl.c

一个关于多重网格的程序

/*  
 * This program is to solve a nonlinear 2-D pde using a one-step linearization 
 *  
 *  It's multigrid, but not full MG. That is, I start at the finest level.
 *  It's written in C, but a fairly transparent version.
 *
 *  the equation I am solving is :
 *     (d/dt)u = Nu+f(x,y,t), where N is a slightly non-linear operator.
 *     N = (\nabla)^2 + function (u)  where \nabla is latex for 
 *          (d^2/dx^2 + d^2/dy^2)
 *     (d/dt)u is obviously a single (time) derivative,
 *     and function(u) is the nonlinear function of u, (no derivatives of u).
 *     the function f(x,y,t) has mostly been chosen to test the code for a known
 *     outcome, but can be any 'nice' function.
 *   The system is 2-dimensional, and has periodic boundary conditions.
 *
 *   N can be any function of u that is continuous on the interval. I have only
 *   rigorously tested for powers of u, exponentials and some trigonometric 
 *   functions. I haven't tried anything else.
 *   
 *   The non-linearity is handled by a one-step linearization at the finest
 *    grid level, and the resulting linear eqn is solved using MG techniques.
 *    See, for example, the book 'Multigrid' by Trottenberg, Oosterlee, Schuller
 *    section 5.3.2 and 5.3.3. The exact eqn solved is explained in the
 *    function mgrid.
 *   There is also a newton relaxation step performed at the finest level,
 *    see for example 'Numerical Methods for Differential Equations' by
 *    Ortega and Poole section 4.4 for a nice introduction.
 *
 *   I've mostly written this as a tutorial code. To that end, I've tried to
 *    explain what's going on throughout the code. The reader is expected to 
 *    have a little familiarity with linear MG and finite difference methods.
 *    And, it's not terribly optimized, as I've tried to keep things transparent
 *    I'll try to provide some hints for a faster code, but profile it & go from
 *    there.
 *
 *   The periodic boundary conditions have a tremendous effect throughout the
 *    code. In particular if you wanted to change to Dirichelet or Neumann
 *    boundary conditions
 *    you would need to change all the loops to (i=1;i<n;i++), assuming the 
 *    boundary was at 0 and n. You'd also have to call the boundary conditions
 *    at every time step, if they are functions of time. I've also assumed 
 *    solution on a square, but that is easily generalizable.
 *
 *    I use the Crank-Nicolson algorithm to integrate in time. It solves the
 *     following equation:
 *      [u(t)-u(t-1)]/tau = (1/2)*[Nu(t) + f(x,y,t) + Nu(t-1) + f(x,y,t-1)]
 *      or:
 *     Nu(t) -(2/tau)u(t)= -[Nu(t-1)+(2/tau) u(t-1)+f(x,y,t-1) + f(x,y,t)]  (**)
 *      where t is the timestep we're now on, and t-1 is the previous time.
 *      Thus (**) is the system we're really solving, and everything on 
 *      the right hand side of the equation is known.
 *      Now, the PDE you solve is  Nu -(2/tau) u = F, where F is the right-hand
 *      side of (**). You will need to hand-tune F, in particular f(x,y,t-1) and
 *      f(x,y,t) for every different equation you run.
 *
 *    
 *    If you don't want to solve a time-dependent PDE, but merely wish to solve
 *      a PDE, ie Nu = f, you'll need to change the routine 'funct' to contain
 *      only f, and remove the while (t<timemax) loop in main.
 *
 *    Debugging. There may well be errors left in the code. Most of the problems
 *      I've had had to do with solving the eqns improperly based on the initial
 *      conditions, and the hand-tuning bit in funct, ie properly calculating
 *      f(x,y,t) based on the initial conditions and the nonlinearity.
 *
 *    This program has mostly been written by Sandra Barsky, with help 
 *      (especially with periodic BC's) from Ralph Giles. My understanding 
 *      of MG techniques is imperfect at best, and there may remain errors 
 *      in this code. If you'd like to point out errors, or know of a way to
 *      improve the accuracy of the solution or make any other comments 
 *      (even if only to let me know you've looked at the program!) please 
 *      send to:    [e-mail]
 */    

//
#define invtau (-2./tau)  //crank-nicholson alog. : need this for linear term
#define tau   0.1  // timestep
#define nrep1 1    // number of pre-smoothing steps
#define nrep2 3    // number of post-smoothing steps
#define nmin 16    // minimum number of grid points
#define nmax 64    // maximum number of grid points
#define timemax 100 // maximum number of time steps

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define pi 4.0*atan(1.0) 
#define xmax pi   //  square I'm solving is [-pi:pi] in both x & y directions
#define xmin -pi  //

#define filesize 64


/* --------------------------------------------------------------------*/
/*  function definitions                                               */
/* --------------------------------------------------------------------*/


void interpolate(double  *uc,double  *u,int nc,int nf);
void restrct(double *uc, double *u,int n);
void restrct2(double *uc, double *u,int n);
void funct(double *f, double *u,int n, int t);
void ic(double *u,int n,int t); /* intial condition, not boundary condition */
double error(double *u,double *uback,int n,int *imx,int *jmx);
void pooofle(double *u, int n, double  *f, double *nlin);
void mgrid(double *u, int n, double *f, int t);
void defectnl(double *d, double *u,int n,double *f);
void newton(double *u, int n, double *d, double *v);
void newtonnl(double *u, int n, double *f);
int index2d(int i,int j, int n);
double errordefect(double *d, int n);


/* --------------------------------------------------------------------*/
/* fn is the nonlinear part of the operator N.
 * df is the derivative of (fn) with respect to u
 */
inline double fn( double u){
  return(u-u*u*u);
}
inline double df( double u){
  return(1.-3.*u*u);
}
/* --------------------------------------------------------------------*/

main()
{
  double *areal; /* areal is the variable I'm solving for. will be called u
		    in the functions that are dependent on this on. I just
		    needed to give it a different name here */
  double *freal; /* freal is the right hand side of (**), ie F */
  double *f,x,y,h,dif,sol;
  int n,nmem,t,i,j;
  FILE *f1 ;
  char filename[filesize];

  n=nmax;
  nmem=(n+1)*(n+1);
  t=1;
  f1=fopen("nlmax","w");/* in this file gets printed the maximum of the 
			   calculated and actual solutions. for comparison */

  areal=malloc(nmem*sizeof(double));
  if(areal==NULL)printf(" oopsies can't allocate areal %d \n",__LINE__);
  ic(areal,n,t);//call initial condition if it's the first time through
	 
  freal=malloc(nmem*sizeof(double));
  if(freal==NULL)printf(" oopsies can't allocate freal %d \n",__LINE__);


  while(t < timemax){

    funct(freal,areal,n,t);
    mgrid(areal,nmax,freal,t);
	
    h=(xmax-xmin)/(double) n;

    fprintf(f1,"%d %lf %lf \n",t,areal[index2d(n/4,n/4,n)],sin(t*tau));
    fflush(f1);

    t++;
  }/* while */
      

  fclose(f1);
  free(areal);
  free(freal);
  
}/* main */

void mgrid(double *u, int n,double *f,int t) 
{
  double *v,*d,*vc,*dc, *uf, *uc,*fc;
  double *uback,*uback2;
  double diffmx,diff,difv,x,y,tol,h;
  double dif,mxd,sumu,sumu2,sum3,sol;
  int i,j,k,imx,jmx,in,nc,irep;
  int ncmem,nsq,nmem;
  char filename[filesize];
  FILE *f1, *f2, *out, *f3, *f8;

  tol=10.e-12; // tol= tolerance. How accurately you want to solve the eqn


  diffmx=tol*5.;
  irep=0;

  nc=n/2;
  nmem=(n+1)*(n+1);  // bigger than [0:n]*[0:n], just to ensure enough space
  ncmem=(nc+1)*(nc+1);// same for the grid that is half-size
  nsq=n*n;  


/* --- requesting space for variables ------------------------------*/

  v=malloc(nmem*sizeof(double));
  if(v==NULL) {printf(" oopsies can't allocate v %d\n",__LINE__);}
  d=malloc(nmem*sizeof(double));
  if(d==NULL) printf(" oopsies can't allocate d %d\n",__LINE__);
  uback=malloc(nmem*sizeof(double));
  if(uback==NULL) printf(" oopsies can't allocate uback %d\n",__LINE__);
  uback2=malloc(nmem*sizeof(double));
  if(uback2==NULL) printf(" oopsies can't allocate uback2  %d\n",__LINE__);
  vc=malloc(ncmem*sizeof(double));
  if(vc==NULL) {printf(" oopsies can't allocate vc  %d\n",__LINE__);}
  uc=malloc(ncmem*sizeof(double));
  if(uc==NULL) {printf(" oopsies can't allocate uc  %d\n",__LINE__);}
  dc=malloc(ncmem*sizeof(double));
  if(dc==NULL) printf(" oopsies can't allocate dc  %d\n",__LINE__);


  for(i=0;i<nmem;i++){  // initialize to 0
    v[i]=0.;
    uback[i]=0.;
    d[i]=0.; }  // end of v,d,uback initialization
  for(i=0;i<ncmem;i++){
    vc[i]=0.;
    dc[i]=0.; }  // end of vc,dc



  h=(xmax-xmin)/(double) n;
  for(i=0;i<n;i++){
    x=xmin+h*i;
    for(j=0;j<n;j++){
      y=xmin+h*j;
      uback2[index2d(i,j,n)]=sin(x)*sin(y)*sin(t*tau);
    }//y     // the actual solution, to compare
  }//x

  diffmx=5*tol; /* tol determines solution tolerance. diffmx will measure 
                   that difference */

  irep=0;  // irep will be the number of MG steps required to solve


  /* +++++++++++++++ start of while loop ++++++++++++++++++*/

  while(diffmx>tol){

    diffmx=0.;   // need to zero this now, will measure later
    irep=irep+1; // counts number of MG steps

		 

    for( in=1;in<=nrep1;in++){  // presmoothing
      newtonnl(u,n,f);
    }//nrep1

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

     /* defectnl calculates the residual. ie, given Nu=f
      *  d=f-Nu, or a measure of how far away you are from the solution */
        
    defectnl(d,u,n,f);

     /* restrict u, have a separate routine for it, in case you want to 
      * use a different restriction for than for other variables */

     restrct2(uc,u,nc);		 


    /* v is the linearlized version of u
     * given, Nu=f
     * define an operator K such that Kv=d, 
     * then u := u+v  (ie new u = old u+ v, v found from above
     * this can be derived from standard Newton's method:
     *  u := u - (Nu -f)/K
     *  where K=(d/du)N
     *
     *  so, if Nu = (\nabla)^2 u + g(u)
     *  K = (\nabla)^2 + (d/du)g
     *
     *  thus, the eqn we now solve is
     *
     *  Kv = (\nabla)^2v + (d/du)g v =d  (#)
     *
     *  and the next iterate of u is u+v
     *
     *  I apply MG techniques to (#)
     */

    restrct(vc,v,nc); 
    restrct(dc,d,nc); // dc is defect on coarse grid

    pooofle(vc,nc,dc,uc); /* funky multigrid part!
			     this solves for vc    */

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

    for(i=1;i<2;i++){
      newton(u,n,d,v);   //   this is 'postsmoothing' for v !!
    }

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

    for( i=0;i<nrep2;i++){
      newtonnl(u,n,f);   	// postsmoothing on u
    }//nrep2


    /* how to tell when you're done?
     * the classic way is to calculate the defect, and when the norm
     * of the defect is less than some predetermined accuracy (tolerance)
     * you're done.
     *
     * I have found that I get better accuracy - ie the solution remains
     * correct farther into time when I use as my condition the difference
     * between the current u, and the previous value of u. If anyone has
     * any ideas for this, please let me know !
     */

    defectnl(d,u,n,f);
    diff=errordefect(d,n)/n/n;

    diffmx=error(u,uback,n,&imx,&jmx); 
    /* diffmx is the condition to test against. for termination of
     * calculation based on defect, do:
     * diffmx=errordefect(d,n)/n/n;
     * diff=error(u,uback,n,&imx,&jmx);
     */ 


    mxd=-10.;
    sum3=0.;
    for(i=0;i<n;i++){
      for(j=0;j<n;j++){
        dif=uback2[index2d(i,j,n)]-u[index2d(i,j,n)]; 
	sum3+=u[index2d(i,j,n)];
        if(dif < 0.e-10) dif= -dif;  
        if(dif > mxd) {mxd= dif;  // maximum difference between actual &
	imx=i;                    // calculated solution. for fun
	jmx=j;                   // imx & jmx calculate where this happens
	} // end if
      }//j
    }//i

     

    if(irep>2000) break; // some break condition
  } /* while */

  printf("%d %d %1.7e %1.7e %1.7e %d %d \n",t,irep,diffmx,diff,mxd,imx,jmx);
  fflush(stdout);
      
      
  /*
   * this will make a file with the x,y coordinates, the calculated &
   * actual solutions & their difference.
   *
   * uncomment if you want this printed out 
   *
  snprintf(filename,filesize,"solution.%d",t);
  f1=fopen(filename,"w");
  for(i=0;i<n;i++){
    x=xmin+h*i;
    for(j=0;j<n;j++){
      y=xmin+h*j;
      sol=sin(x)*sin(y)*sin(t*tau);
      dif=sol-u[index2d(i,j,n)];
      fprintf(f1,"%lf %lf %lf %lf %lf\n",x,y,u[index2d(i,j,n)],sol,dif,uback[index2d(i,j,n)]);
    }
  }
  fclose(f1);

  */

  free(uback);
  free(uback2);