Actual source code: ex3.c
petsc-3.4.5 2014-06-29
2: static char help[] = "Solves 1D heat equation with FEM formulation.\n\
3: Input arguments are\n\
4: -useAlhs: solve Alhs*U' = (Arhs*U + g) \n\
5: otherwise, solve U' = inv(Alhs)*(Arhs*U + g) \n\n";
7: /*--------------------------------------------------------------------------
8: Solves 1D heat equation U_t = U_xx with FEM formulation:
9: Alhs*U' = rhs (= Arhs*U + g)
10: We thank Chris Cox <clcox@clemson.edu> for contributing the original code
11: ----------------------------------------------------------------------------*/
13: #include <petscksp.h>
14: #include <petscts.h>
16: /* special variable - max size of all arrays */
17: #define num_z 60
19: /*
20: User-defined application context - contains data needed by the
21: application-provided call-back routines.
22: */
23: typedef struct {
24: Mat Amat; /* left hand side matrix */
25: Vec ksp_rhs,ksp_sol; /* working vectors for formulating inv(Alhs)*(Arhs*U+g) */
26: int max_probsz; /* max size of the problem */
27: PetscBool useAlhs; /* flag (1 indicates solving Alhs*U' = Arhs*U+g */
28: int nz; /* total number of grid points */
29: PetscInt m; /* total number of interio grid points */
30: Vec solution; /* global exact ts solution vector */
31: PetscScalar *z; /* array of grid points */
32: PetscBool debug; /* flag (1 indicates activation of debugging printouts) */
33: } AppCtx;
35: extern PetscScalar exact(PetscScalar,PetscReal);
36: extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);
37: extern PetscErrorCode Petsc_KSPSolve(AppCtx*);
38: extern PetscScalar bspl(PetscScalar*,PetscScalar,PetscInt,PetscInt,PetscInt[][2],PetscInt);
39: extern void femBg(PetscScalar[][3],PetscScalar*,PetscInt,PetscScalar*,PetscReal);
40: extern void femA(AppCtx*,PetscInt,PetscScalar*);
41: extern void rhs(AppCtx*,PetscScalar*, PetscInt, PetscScalar*,PetscReal);
42: extern PetscErrorCode RHSfunction(TS,PetscReal,Vec,Vec,void*);
46: int main(int argc,char **argv)
47: {
48: PetscInt i,m,nz,steps,max_steps,k,nphase=1;
49: PetscScalar zInitial,zFinal,val,*z;
50: PetscReal stepsz[4],T,ftime;
52: TS ts;
53: SNES snes;
54: Mat Jmat;
55: AppCtx appctx; /* user-defined application context */
56: Vec init_sol; /* ts solution vector */
57: PetscMPIInt size;
59: PetscInitialize(&argc,&argv,(char*)0,help);
60: MPI_Comm_size(PETSC_COMM_WORLD,&size);
61: if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only");
63: PetscOptionsGetInt(NULL,"-nphase",&nphase,NULL);
64: if (nphase > 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nphase must be an integer between 1 and 3");
66: /* initializations */
67: zInitial = 0.0;
68: zFinal = 1.0;
69: T = 0.014/nphase;
70: nz = num_z;
71: m = nz-2;
72: appctx.nz = nz;
73: max_steps = (PetscInt)10000;
75: appctx.m = m;
76: appctx.max_probsz = nz;
77: appctx.debug = PETSC_FALSE;
78: appctx.useAlhs = PETSC_FALSE;
80: PetscOptionsHasName(NULL,"-debug",&appctx.debug);
81: PetscOptionsHasName(NULL,"-useAlhs",&appctx.useAlhs);
83: /* create vector to hold ts solution */
84: /*-----------------------------------*/
85: VecCreate(PETSC_COMM_WORLD, &init_sol);
86: VecSetSizes(init_sol, PETSC_DECIDE, m);
87: VecSetFromOptions(init_sol);
89: /* create vector to hold true ts soln for comparison */
90: VecDuplicate(init_sol, &appctx.solution);
92: /* create LHS matrix Amat */
93: /*------------------------*/
94: MatCreateSeqAIJ(PETSC_COMM_WORLD, m, m, 3, NULL, &appctx.Amat);
95: MatSetFromOptions(appctx.Amat);
96: MatSetUp(appctx.Amat);
97: /* set space grid points - interio points only! */
98: PetscMalloc((nz+1)*sizeof(PetscScalar),&z);
99: for (i=0; i<nz; i++) z[i]=(i)*((zFinal-zInitial)/(nz-1));
100: appctx.z = z;
101: femA(&appctx,nz,z);
103: /* create the jacobian matrix */
104: /*----------------------------*/
105: MatCreate(PETSC_COMM_WORLD, &Jmat);
106: MatSetSizes(Jmat,PETSC_DECIDE,PETSC_DECIDE,m,m);
107: MatSetFromOptions(Jmat);
108: MatSetUp(Jmat);
110: /* create working vectors for formulating rhs=inv(Alhs)*(Arhs*U + g) */
111: VecDuplicate(init_sol,&appctx.ksp_rhs);
112: VecDuplicate(init_sol,&appctx.ksp_sol);
114: /* set intial guess */
115: /*------------------*/
116: for (i=0; i<nz-2; i++) {
117: val = exact(z[i+1], 0.0);
118: VecSetValue(init_sol,i,(PetscScalar)val,INSERT_VALUES);
119: }
120: VecAssemblyBegin(init_sol);
121: VecAssemblyEnd(init_sol);
123: /*create a time-stepping context and set the problem type */
124: /*--------------------------------------------------------*/
125: TSCreate(PETSC_COMM_WORLD, &ts);
126: TSSetProblemType(ts,TS_NONLINEAR);
128: /* set time-step method */
129: TSSetType(ts,TSCN);
131: /* Set optional user-defined monitoring routine */
132: TSMonitorSet(ts,Monitor,&appctx,NULL);
133: /* set the right hand side of U_t = RHSfunction(U,t) */
134: TSSetRHSFunction(ts,NULL,(PetscErrorCode (*)(TS,PetscScalar,Vec,Vec,void*))RHSfunction,&appctx);
136: if (appctx.useAlhs) {
137: /* set the left hand side matrix of Amat*U_t = rhs(U,t) */
139: /* Note: this approach is incompatible with the finite differenced Jacobian set below because we can't restore the
140: * Alhs matrix without making a copy. Either finite difference the entire thing or use analytic Jacobians in both
141: * places.
142: */
143: TSSetIFunction(ts,NULL,TSComputeIFunctionLinear,&appctx);
144: TSSetIJacobian(ts,appctx.Amat,appctx.Amat,TSComputeIJacobianConstant,&appctx);
145: }
147: /* use petsc to compute the jacobian by finite differences */
148: TSGetSNES(ts,&snes);
149: SNESSetJacobian(snes,Jmat,Jmat,SNESComputeJacobianDefault,NULL);
151: /* get the command line options if there are any and set them */
152: TSSetFromOptions(ts);
154: #if defined(PETSC_HAVE_SUNDIALS)
155: {
156: TSType type;
157: PetscBool sundialstype=PETSC_FALSE;
158: TSGetType(ts,&type);
159: PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&sundialstype);
160: if (sundialstype && appctx.useAlhs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use Alhs formulation for TSSUNDIALS type");
161: }
162: #endif
163: /* Sets the initial solution */
164: TSSetSolution(ts,init_sol);
166: stepsz[0] = 1.0/(2.0*(nz-1)*(nz-1)); /* (mesh_size)^2/2.0 */
167: ftime = 0.0;
168: for (k=0; k<nphase; k++) {
169: if (nphase > 1) printf("Phase %d: initial time %g, stepsz %g, duration: %g\n",k,ftime,stepsz[k],(k+1)*T);
170: TSSetInitialTimeStep(ts,ftime,stepsz[k]);
171: TSSetDuration(ts,max_steps,(k+1)*T);
173: /* loop over time steps */
174: /*----------------------*/
175: TSSolve(ts,init_sol);
176: TSGetSolveTime(ts,&ftime);
177: TSGetTimeStepNumber(ts,&steps);
178: stepsz[k+1] = stepsz[k]*1.5; /* change step size for the next phase */
179: }
181: /* free space */
182: TSDestroy(&ts);
183: MatDestroy(&appctx.Amat);
184: MatDestroy(&Jmat);
185: VecDestroy(&appctx.ksp_rhs);
186: VecDestroy(&appctx.ksp_sol);
187: VecDestroy(&init_sol);
188: VecDestroy(&appctx.solution);
189: PetscFree(z);
191: PetscFinalize();
192: return 0;
193: }
195: /*------------------------------------------------------------------------
196: Set exact solution
197: u(z,t) = sin(6*PI*z)*exp(-36.*PI*PI*t) + 3.*sin(2*PI*z)*exp(-4.*PI*PI*t)
198: --------------------------------------------------------------------------*/
199: PetscScalar exact(PetscScalar z,PetscReal t)
200: {
201: PetscScalar val, ex1, ex2;
203: ex1 = exp(-36.*PETSC_PI*PETSC_PI*t);
204: ex2 = exp(-4.*PETSC_PI*PETSC_PI*t);
205: val = sin(6*PETSC_PI*z)*ex1 + 3.*sin(2*PETSC_PI*z)*ex2;
206: return val;
207: }
211: /*
212: Monitor - User-provided routine to monitor the solution computed at
213: each timestep. This example plots the solution and computes the
214: error in two different norms.
216: Input Parameters:
217: ts - the timestep context
218: step - the count of the current step (with 0 meaning the
219: initial condition)
220: time - the current time
221: u - the solution at this timestep
222: ctx - the user-provided context for this monitoring routine.
223: In this case we use the application context which contains
224: information about the problem size, workspace and the exact
225: solution.
226: */
227: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec u,void *ctx)
228: {
229: AppCtx *appctx = (AppCtx*)ctx;
231: PetscInt i,m=appctx->m;
232: PetscReal norm_2,norm_max,h=1.0/(m+1);
233: PetscScalar *u_exact;
235: /* Compute the exact solution */
236: VecGetArray(appctx->solution,&u_exact);
237: for (i=0; i<m; i++) u_exact[i] = exact(appctx->z[i+1],time);
238: VecRestoreArray(appctx->solution,&u_exact);
240: /* Print debugging information if desired */
241: if (appctx->debug) {
242: PetscPrintf(PETSC_COMM_SELF,"Computed solution vector at time %g\n",time);
243: VecView(u,PETSC_VIEWER_STDOUT_SELF);
244: PetscPrintf(PETSC_COMM_SELF,"Exact solution vector\n");
245: VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF);
246: }
248: /* Compute the 2-norm and max-norm of the error */
249: VecAXPY(appctx->solution,-1.0,u);
250: VecNorm(appctx->solution,NORM_2,&norm_2);
252: norm_2 = PetscSqrtReal(h)*norm_2;
253: VecNorm(appctx->solution,NORM_MAX,&norm_max);
255: PetscPrintf(PETSC_COMM_SELF,"Timestep %D: time = %G, 2-norm error = %6.4f, max norm error = %6.4f\n",
256: step,time,norm_2,norm_max);
258: /*
259: Print debugging information if desired
260: */
261: if (appctx->debug) {
262: PetscPrintf(PETSC_COMM_SELF,"Error vector\n");
263: VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF);
264: }
265: return 0;
266: }
268: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
269: %% Function to solve a linear system using KSP %%
270: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
272: PetscErrorCode Petsc_KSPSolve(AppCtx *obj)
273: {
275: KSP ksp;
276: PC pc;
278: /*create the ksp context and set the operators,that is, associate the system matrix with it*/
279: KSPCreate(PETSC_COMM_WORLD,&ksp);
280: KSPSetOperators(ksp,obj->Amat,obj->Amat,DIFFERENT_NONZERO_PATTERN);
282: /*get the preconditioner context, set its type and the tolerances*/
283: KSPGetPC(ksp,&pc);
284: PCSetType(pc,PCLU);
285: KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
287: /*get the command line options if there are any and set them*/
288: KSPSetFromOptions(ksp);
290: /*get the linear system (ksp) solve*/
291: KSPSolve(ksp,obj->ksp_rhs,obj->ksp_sol);
293: KSPDestroy(&ksp);
294: return 0;
295: }
297: /***********************************************************************
298: * Function to return value of basis function or derivative of basis *
299: * function. *
300: ***********************************************************************
301: * *
302: * Arguments: *
303: * x = array of xpoints or nodal values *
304: * xx = point at which the basis function is to be *
305: * evaluated. *
306: * il = interval containing xx. *
307: * iq = indicates which of the two basis functions in *
308: * interval intrvl should be used *
309: * nll = array containing the endpoints of each interval. *
310: * id = If id ~= 2, the value of the basis function *
311: * is calculated; if id = 2, the value of the *
312: * derivative of the basis function is returned. *
313: ***********************************************************************/
315: PetscScalar bspl(PetscScalar *x, PetscScalar xx,PetscInt il,PetscInt iq,PetscInt nll[][2],PetscInt id)
316: {
317: PetscScalar x1,x2,bfcn;
318: PetscInt i1,i2,iq1,iq2;
320: /*** Determine which basis function in interval intrvl is to be used in ***/
321: iq1 = iq;
322: if (iq1==0) iq2 = 1;
323: else iq2 = 0;
325: /*** Determine endpoint of the interval intrvl ***/
326: i1=nll[il][iq1];
327: i2=nll[il][iq2];
329: /*** Determine nodal values at the endpoints of the interval intrvl ***/
330: x1=x[i1];
331: x2=x[i2];
332: /* printf("x1=%g\tx2=%g\txx=%g\n",x1,x2,xx); */
333: /*** Evaluate basis function ***/
334: if (id == 2) bfcn=(1.0)/(x1-x2);
335: else bfcn=(xx-x2)/(x1-x2);
336: /* printf("bfcn=%g\n",bfcn); */
337: return bfcn;
338: }
340: /*---------------------------------------------------------
341: Function called by rhs function to get B and g
342: ---------------------------------------------------------*/
343: void femBg(PetscScalar btri[][3],PetscScalar *f,PetscInt nz,PetscScalar *z, PetscReal t)
344: {
345: PetscInt i,j,jj,il,ip,ipp,ipq,iq,iquad,iqq;
346: PetscInt nli[num_z][2],indx[num_z];
347: PetscScalar dd,dl,zip,zipq,zz,bb,b_z,bbb,bb_z,bij;
348: PetscScalar zquad[num_z][3],dlen[num_z],qdwt[3];
350: /* initializing everything - btri and f are initialized in rhs.c */
351: for (i=0; i < nz; i++) {
352: nli[i][0] = 0;
353: nli[i][1] = 0;
354: indx[i] = 0;
355: zquad[i][0] = 0.0;
356: zquad[i][1] = 0.0;
357: zquad[i][2] = 0.0;
358: dlen[i] = 0.0;
359: } /*end for (i)*/
361: /* quadrature weights */
362: qdwt[0] = 1.0/6.0;
363: qdwt[1] = 4.0/6.0;
364: qdwt[2] = 1.0/6.0;
366: /* 1st and last nodes have Dirichlet boundary condition -
367: set indices there to -1 */
369: for (i=0; i < nz-1; i++) indx[i] = i-1;
370: indx[nz-1] = -1;
372: ipq = 0;
373: for (il=0; il < nz-1; il++) {
374: ip = ipq;
375: ipq = ip+1;
376: zip = z[ip];
377: zipq = z[ipq];
378: dl = zipq-zip;
379: zquad[il][0] = zip;
380: zquad[il][1] = (0.5)*(zip+zipq);
381: zquad[il][2] = zipq;
382: dlen[il] = fabs(dl);
383: nli[il][0] = ip;
384: nli[il][1] = ipq;
385: }
387: for (il=0; il < nz-1; il++) {
388: for (iquad=0; iquad < 3; iquad++) {
389: dd = (dlen[il])*(qdwt[iquad]);
390: zz = zquad[il][iquad];
392: for (iq=0; iq < 2; iq++) {
393: ip = nli[il][iq];
394: bb = bspl(z,zz,il,iq,nli,1);
395: b_z = bspl(z,zz,il,iq,nli,2);
396: i = indx[ip];
398: if (i > -1) {
399: for (iqq=0; iqq < 2; iqq++) {
400: ipp = nli[il][iqq];
401: bbb = bspl(z,zz,il,iqq,nli,1);
402: bb_z = bspl(z,zz,il,iqq,nli,2);
403: j = indx[ipp];
404: bij = -b_z*bb_z;
406: if (j > -1) {
407: jj = 1+j-i;
408: btri[i][jj] += bij*dd;
409: } else {
410: f[i] += bij*dd*exact(z[ipp], t);
411: /* f[i] += 0.0; */
412: /* if (il==0 && j==-1) { */
413: /* f[i] += bij*dd*exact(zz,t); */
414: /* }*/ /*end if*/
415: } /*end else*/
416: } /*end for (iqq)*/
417: } /*end if (i>0)*/
418: } /*end for (iq)*/
419: } /*end for (iquad)*/
420: } /*end for (il)*/
421: return;
422: }
424: void femA(AppCtx *obj,PetscInt nz,PetscScalar *z)
425: {
426: PetscInt i,j,il,ip,ipp,ipq,iq,iquad,iqq;
427: PetscInt nli[num_z][2],indx[num_z];
428: PetscScalar dd,dl,zip,zipq,zz,bb,bbb,aij;
429: PetscScalar rquad[num_z][3],dlen[num_z],qdwt[3],add_term;
432: /* initializing everything */
434: for (i=0; i < nz; i++) {
435: nli[i][0] = 0;
436: nli[i][1] = 0;
437: indx[i] = 0;
438: rquad[i][0] = 0.0;
439: rquad[i][1] = 0.0;
440: rquad[i][2] = 0.0;
441: dlen[i] = 0.0;
442: } /*end for (i)*/
444: /* quadrature weights */
445: qdwt[0] = 1.0/6.0;
446: qdwt[1] = 4.0/6.0;
447: qdwt[2] = 1.0/6.0;
449: /* 1st and last nodes have Dirichlet boundary condition -
450: set indices there to -1 */
452: for (i=0; i < nz-1; i++) indx[i]=i-1;
453: indx[nz-1]=-1;
455: ipq = 0;
457: for (il=0; il < nz-1; il++) {
458: ip = ipq;
459: ipq = ip+1;
460: zip = z[ip];
461: zipq = z[ipq];
462: dl = zipq-zip;
463: rquad[il][0] = zip;
464: rquad[il][1] = (0.5)*(zip+zipq);
465: rquad[il][2] = zipq;
466: dlen[il] = fabs(dl);
467: nli[il][0] = ip;
468: nli[il][1] = ipq;
469: } /*end for (il)*/
471: for (il=0; il < nz-1; il++) {
472: for (iquad=0; iquad < 3; iquad++) {
473: dd = (dlen[il])*(qdwt[iquad]);
474: zz = rquad[il][iquad];
476: for (iq=0; iq < 2; iq++) {
477: ip = nli[il][iq];
478: bb = bspl(z,zz,il,iq,nli,1);
479: i = indx[ip];
480: if (i > -1) {
481: for (iqq=0; iqq < 2; iqq++) {
482: ipp = nli[il][iqq];
483: bbb = bspl(z,zz,il,iqq,nli,1);
484: j = indx[ipp];
485: aij = bb*bbb;
486: if (j > -1) {
487: add_term = aij*dd;
488: MatSetValue(obj->Amat,i,j,add_term,ADD_VALUES);
489: }/*endif*/
490: } /*end for (iqq)*/
491: } /*end if (i>0)*/
492: } /*end for (iq)*/
493: } /*end for (iquad)*/
494: } /*end for (il)*/
495: MatAssemblyBegin(obj->Amat,MAT_FINAL_ASSEMBLY);
496: MatAssemblyEnd(obj->Amat,MAT_FINAL_ASSEMBLY);
497: return;
498: }
500: /*---------------------------------------------------------
501: Function to fill the rhs vector with
502: By + g values ****
503: ---------------------------------------------------------*/
504: void rhs(AppCtx *obj,PetscScalar *y, PetscInt nz, PetscScalar *z, PetscReal t)
505: {
506: PetscInt i,j,js,je,jj;
507: PetscScalar val,g[num_z],btri[num_z][3],add_term;
510: for (i=0; i < nz-2; i++) {
511: for (j=0; j <= 2; j++) btri[i][j]=0.0;
512: g[i] = 0.0;
513: }
515: /* call femBg to set the tri-diagonal b matrix and vector g */
516: femBg(btri,g,nz,z,t);
518: /* setting the entries of the right hand side vector */
519: for (i=0; i < nz-2; i++) {
520: val = 0.0;
521: js = 0;
522: if (i == 0) js = 1;
523: je = 2;
524: if (i == nz-2) je = 1;
526: for (jj=js; jj <= je; jj++) {
527: j = i+jj-1;
528: val += (btri[i][jj])*(y[j]);
529: }
530: add_term = val + g[i];
531: VecSetValue(obj->ksp_rhs,(PetscInt)i,(PetscScalar)add_term,INSERT_VALUES);
532: }
533: VecAssemblyBegin(obj->ksp_rhs);
534: VecAssemblyEnd(obj->ksp_rhs);
536: /* return to main driver function */
537: return;
538: }
540: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
541: %% Function to form the right hand side of the time-stepping problem. %%
542: %% -------------------------------------------------------------------------------------------%%
543: if (useAlhs):
544: globalout = By+g
545: else if (!useAlhs):
546: globalout = f(y,t)=Ainv(By+g),
547: in which the ksp solver to transform the problem A*ydot=By+g
548: to the problem ydot=f(y,t)=inv(A)*(By+g)
549: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
551: PetscErrorCode RHSfunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
552: {
554: AppCtx *obj = (AppCtx*)ctx;
555: PetscScalar *soln_ptr,soln[num_z-2];
556: PetscInt i,nz=obj->nz;
557: PetscReal time;
559: /* get the previous solution to compute updated system */
560: VecGetArray(globalin,&soln_ptr);
561: for (i=0; i < num_z-2; i++) soln[i] = soln_ptr[i];
562: VecRestoreArray(globalin,&soln_ptr);
564: /* clear out the matrix and rhs for ksp to keep things straight */
565: VecSet(obj->ksp_rhs,(PetscScalar)0.0);
567: time = t;
568: /* get the updated system */
569: rhs(obj,soln,nz,obj->z,time); /* setup of the By+g rhs */
571: /* do a ksp solve to get the rhs for the ts problem */
572: if (obj->useAlhs) {
573: /* ksp_sol = ksp_rhs */
574: VecCopy(obj->ksp_rhs,globalout);
575: } else {
576: /* ksp_sol = inv(Amat)*ksp_rhs */
577: Petsc_KSPSolve(obj);
578: VecCopy(obj->ksp_sol,globalout);
579: }
580: return 0;
581: }