/* * 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 #include #include #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;itol){ 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 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 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 nmin){ for(i=0;i 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=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 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