Actual source code: ex3.c
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;
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: Mat Jmat;
54: AppCtx appctx; /* user-defined application context */
55: Vec init_sol; /* ts solution vector */
56: PetscMPIInt size;
58: PetscInitialize(&argc,&argv,(char*)0,help);
59: MPI_Comm_size(PETSC_COMM_WORLD,&size);
60: if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only");
62: PetscOptionsHasName(PETSC_NULL,"-debug",&appctx.debug);
63: PetscOptionsHasName(PETSC_NULL,"-useAlhs",&appctx.useAlhs);
64: PetscOptionsGetInt(PETSC_NULL,"-nphase",&nphase,PETSC_NULL);
65: if (nphase > 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nphase must be an integer between 1 and 3");
67: /* initializations */
68: zInitial = 0.0;
69: zFinal = 1.0;
70: T = 0.014/nphase;
71: nz = num_z;
72: m = nz-2;
73: appctx.nz = nz;
74: max_steps = (PetscInt)10000;
76: appctx.m = m;
77: appctx.max_probsz = nz;
78: appctx.debug = PETSC_FALSE;
79: appctx.useAlhs = PETSC_FALSE;
81: /* create vector to hold ts solution */
82: /*-----------------------------------*/
83: VecCreate(PETSC_COMM_WORLD, &init_sol);
84: VecSetSizes(init_sol, PETSC_DECIDE, m);
85: VecSetFromOptions(init_sol);
87: /* create vector to hold true ts soln for comparison */
88: VecDuplicate(init_sol, &appctx.solution);
90: /* create LHS matrix Amat */
91: /*------------------------*/
92: MatCreateSeqAIJ(PETSC_COMM_WORLD, m, m, 3, PETSC_NULL, &appctx.Amat);
93: MatSetFromOptions(appctx.Amat);
94: /* set space grid points - interio points only! */
95: PetscMalloc((nz+1)*sizeof(PetscScalar),&z);
96: for (i=0; i<nz; i++) z[i]=(i)*((zFinal-zInitial)/(nz-1));
97: appctx.z = z;
98: femA(&appctx,nz,z);
100: /* create the jacobian matrix */
101: /*----------------------------*/
102: MatCreate(PETSC_COMM_WORLD, &Jmat);
103: MatSetSizes(Jmat,PETSC_DECIDE,PETSC_DECIDE,m,m);
104: MatSetFromOptions(Jmat);
106: /* create working vectors for formulating rhs=inv(Alhs)*(Arhs*U + g) */
107: VecDuplicate(init_sol,&appctx.ksp_rhs);
108: VecDuplicate(init_sol,&appctx.ksp_sol);
110: /* set intial guess */
111: /*------------------*/
112: for(i=0; i<nz-2; i++){
113: val = exact(z[i+1], 0.0);
114: VecSetValue(init_sol,i,(PetscScalar)val,INSERT_VALUES);
115: }
116: VecAssemblyBegin(init_sol);
117: VecAssemblyEnd(init_sol);
119: /*create a time-stepping context and set the problem type */
120: /*--------------------------------------------------------*/
121: TSCreate(PETSC_COMM_WORLD, &ts);
122: TSSetProblemType(ts,TS_NONLINEAR);
124: /* set time-step method */
125: TSSetType(ts,TSCN);
127: /* Set optional user-defined monitoring routine */
128: TSMonitorSet(ts,Monitor,&appctx,PETSC_NULL);
129: /* set the right hand side of U_t = RHSfunction(U,t) */
130: TSSetRHSFunction(ts,PETSC_NULL,(PetscErrorCode (*)(TS,PetscScalar,Vec,Vec,void*))RHSfunction,&appctx);
132: if (appctx.useAlhs){
133: /* set the left hand side matrix of Amat*U_t = rhs(U,t) */
134: TSSetIFunction(ts,PETSC_NULL,TSComputeIFunctionLinear,&appctx);
135: TSSetIJacobian(ts,appctx.Amat,appctx.Amat,TSComputeIJacobianConstant,&appctx);
136: }
138: /* use petsc to compute the jacobian by finite differences */
139: TSSetRHSJacobian(ts,Jmat,Jmat,TSDefaultComputeJacobian,&appctx);
141: /* get the command line options if there are any and set them */
142: TSSetFromOptions(ts);
144: #ifdef PETSC_HAVE_SUNDIALS
145: {
146: const TSType type;
147: PetscBool sundialstype=PETSC_FALSE;
148: TSGetType(ts,&type);
149: PetscTypeCompare((PetscObject)ts,TSSUNDIALS,&sundialstype);
150: if (sundialstype && appctx.useAlhs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use Alhs formulation for TSSUNDIALS type");
151: }
152: #endif
153: /* Sets the initial solution */
154: TSSetSolution(ts,init_sol);
156: stepsz[0] = 1.0/(2.0*(nz-1)*(nz-1)); /* (mesh_size)^2/2.0 */
157: ftime = 0.0;
158: for (k=0; k<nphase; k++){
159: if (nphase > 1) {
160: printf("Phase %d: initial time %g, stepsz %g, duration: %g\n",k,ftime,stepsz[k],(k+1)*T);
161: }
162: TSSetInitialTimeStep(ts,ftime,stepsz[k]);
163: TSSetDuration(ts,max_steps,(k+1)*T);
165: /* loop over time steps */
166: /*----------------------*/
167: TSSolve(ts,init_sol,&ftime);
168: TSGetTimeStepNumber(ts,&steps);
169: stepsz[k+1] = stepsz[k]*1.5; /* change step size for the next phase */
170: }
172: /* free space */
173: TSDestroy(&ts);
174: MatDestroy(&appctx.Amat);
175: MatDestroy(&Jmat);
176: VecDestroy(&appctx.ksp_rhs);
177: VecDestroy(&appctx.ksp_sol);
178: VecDestroy(&init_sol);
179: VecDestroy(&appctx.solution);
180: PetscFree(z);
182: PetscFinalize();
183: return 0;
184: }
186: /*------------------------------------------------------------------------
187: Set exact solution
188: u(z,t) = sin(6*PI*z)*exp(-36.*PI*PI*t) + 3.*sin(2*PI*z)*exp(-4.*PI*PI*t)
189: --------------------------------------------------------------------------*/
190: PetscScalar exact(PetscScalar z,PetscReal t)
191: {
192: PetscScalar val, ex1, ex2;
194: ex1 = exp(-36.*PETSC_PI*PETSC_PI*t);
195: ex2 = exp(-4.*PETSC_PI*PETSC_PI*t);
196: val = sin(6*PETSC_PI*z)*ex1 + 3.*sin(2*PETSC_PI*z)*ex2;
197: return val;
198: }
202: /*
203: Monitor - User-provided routine to monitor the solution computed at
204: each timestep. This example plots the solution and computes the
205: error in two different norms.
207: Input Parameters:
208: ts - the timestep context
209: step - the count of the current step (with 0 meaning the
210: initial condition)
211: time - the current time
212: u - the solution at this timestep
213: ctx - the user-provided context for this monitoring routine.
214: In this case we use the application context which contains
215: information about the problem size, workspace and the exact
216: solution.
217: */
218: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec u,void *ctx)
219: {
220: AppCtx *appctx = (AppCtx*)ctx;
222: PetscInt i,m=appctx->m;
223: PetscReal norm_2,norm_max,h=1.0/(m+1);
224: PetscScalar *u_exact;
226: /* Compute the exact solution */
227: VecGetArray(appctx->solution,&u_exact);
228: for (i=0; i<m; i++){
229: u_exact[i] = exact(appctx->z[i+1],time);
230: }
231: VecRestoreArray(appctx->solution,&u_exact);
233: /* Print debugging information if desired */
234: if (appctx->debug) {
235: PetscPrintf(PETSC_COMM_SELF,"Computed solution vector at time %g\n",time);
236: VecView(u,PETSC_VIEWER_STDOUT_SELF);
237: PetscPrintf(PETSC_COMM_SELF,"Exact solution vector\n");
238: VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF);
239: }
241: /* Compute the 2-norm and max-norm of the error */
242: VecAXPY(appctx->solution,-1.0,u);
243: VecNorm(appctx->solution,NORM_2,&norm_2);
245: norm_2 = sqrt(h)*norm_2;
246: VecNorm(appctx->solution,NORM_MAX,&norm_max);
248: PetscPrintf(PETSC_COMM_SELF,"Timestep %D: time = %G, 2-norm error = %6.4f, max norm error = %6.4f\n",
249: step,time,norm_2,norm_max);
251: /*
252: Print debugging information if desired
253: */
254: if (appctx->debug) {
255: PetscPrintf(PETSC_COMM_SELF,"Error vector\n");
256: VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF);
257: }
258: return 0;
259: }
261: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
262: %% Function to solve a linear system using KSP %%
263: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
265: void Petsc_KSPSolve(AppCtx *obj)
266: {
267: PetscErrorCode ierr;
268: KSP ksp;
269: PC pc;
271: /*create the ksp context and set the operators,that is, associate the system matrix with it*/
272: KSPCreate(PETSC_COMM_WORLD,&ksp);
273: KSPSetOperators(ksp,obj->Amat,obj->Amat,DIFFERENT_NONZERO_PATTERN);
275: /*get the preconditioner context, set its type and the tolerances*/
276: KSPGetPC(ksp,&pc);
277: PCSetType(pc,PCLU);
278: KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
280: /*get the command line options if there are any and set them*/
281: KSPSetFromOptions(ksp);
283: /*get the linear system (ksp) solve*/
284: KSPSolve(ksp,obj->ksp_rhs,obj->ksp_sol);
286: KSPDestroy(&ksp);
287: return;
288: }
290: /***********************************************************************
291: * Function to return value of basis function or derivative of basis *
292: * function. *
293: ***********************************************************************
294: * *
295: * Arguments: *
296: * x = array of xpoints or nodal values *
297: * xx = point at which the basis function is to be *
298: * evaluated. *
299: * il = interval containing xx. *
300: * iq = indicates which of the two basis functions in *
301: * interval intrvl should be used *
302: * nll = array containing the endpoints of each interval. *
303: * id = If id ~= 2, the value of the basis function *
304: * is calculated; if id = 2, the value of the *
305: * derivative of the basis function is returned. *
306: ***********************************************************************/
308: PetscScalar bspl(PetscScalar *x, PetscScalar xx,PetscInt il,PetscInt iq,PetscInt nll[][2],PetscInt id)
309: {
310: PetscScalar x1,x2,bfcn;
311: PetscInt i1,i2,iq1,iq2;
313: /*** Determine which basis function in interval intrvl is to be used in ***/
314: iq1 = iq;
315: if(iq1==0) iq2 = 1;
316: else iq2 = 0;
318: /*** Determine endpoint of the interval intrvl ***/
319: i1=nll[il][iq1];
320: i2=nll[il][iq2];
322: /*** Determine nodal values at the endpoints of the interval intrvl ***/
323: x1=x[i1];
324: x2=x[i2];
325: //printf("x1=%g\tx2=%g\txx=%g\n",x1,x2,xx);
326: /*** Evaluate basis function ***/
327: if(id == 2) bfcn=(1.0)/(x1-x2);
328: else bfcn=(xx-x2)/(x1-x2);
329: //printf("bfcn=%g\n",bfcn);
330: return bfcn;
331: }
333: /*---------------------------------------------------------
334: Function called by rhs function to get B and g
335: ---------------------------------------------------------*/
336: void femBg(PetscScalar btri[][3],PetscScalar *f,PetscInt nz,PetscScalar *z, PetscReal t)
337: {
338: PetscInt i,j,jj,il,ip,ipp,ipq,iq,iquad,iqq;
339: PetscInt nli[num_z][2],indx[num_z];
340: PetscScalar dd,dl,zip,zipq,zz,bb,b_z,bbb,bb_z,bij;
341: PetscScalar zquad[num_z][3],dlen[num_z],qdwt[3];
343: /* initializing everything - btri and f are initialized in rhs.c */
344: for(i=0; i < nz; i++){
345: nli[i][0] = 0;
346: nli[i][1] = 0;
347: indx[i] = 0;
348: zquad[i][0] = 0.0;
349: zquad[i][1] = 0.0;
350: zquad[i][2] = 0.0;
351: dlen[i] = 0.0;
352: }/*end for(i)*/
354: /* quadrature weights */
355: qdwt[0] = 1.0/6.0;
356: qdwt[1] = 4.0/6.0;
357: qdwt[2] = 1.0/6.0;
359: /* 1st and last nodes have Dirichlet boundary condition -
360: set indices there to -1 */
362: for(i=0; i < nz-1; i++){
363: indx[i]=i-1;
364: }
365: indx[nz-1]=-1;
367: ipq = 0;
368: for (il=0; il < nz-1; il++){
369: ip = ipq;
370: ipq = ip+1;
371: zip = z[ip];
372: zipq = z[ipq];
373: dl = zipq-zip;
374: zquad[il][0] = zip;
375: zquad[il][1] = (0.5)*(zip+zipq);
376: zquad[il][2] = zipq;
377: dlen[il] = fabs(dl);
378: nli[il][0] = ip;
379: nli[il][1] = ipq;
380: }
382: for (il=0; il < nz-1; il++){
383: for (iquad=0; iquad < 3; iquad++){
384: dd = (dlen[il])*(qdwt[iquad]);
385: zz = zquad[il][iquad];
387: for (iq=0; iq < 2; iq++){
388: ip = nli[il][iq];
389: bb = bspl(z,zz,il,iq,nli,1);
390: b_z = bspl(z,zz,il,iq,nli,2);
391: i = indx[ip];
393: if(i > -1){
394: for(iqq=0; iqq < 2; iqq++){
395: ipp = nli[il][iqq];
396: bbb = bspl(z,zz,il,iqq,nli,1);
397: bb_z = bspl(z,zz,il,iqq,nli,2);
398: j = indx[ipp];
399: bij = -b_z*bb_z;
401: if (j > -1){
402: jj = 1+j-i;
403: btri[i][jj] += bij*dd;
404: } else {
405: f[i] += bij*dd*exact(z[ipp], t);
406: // f[i] += 0.0;
407: // if(il==0 && j==-1){
408: // f[i] += bij*dd*exact(zz,t);
409: // }/*end if*/
410: } /*end else*/
411: }/*end for(iqq)*/
412: }/*end if(i>0)*/
413: }/*end for(iq)*/
414: }/*end for(iquad)*/
415: }/*end for(il)*/
416: return;
417: }
419: void femA(AppCtx *obj,PetscInt nz,PetscScalar *z)
420: {
421: PetscInt i,j,il,ip,ipp,ipq,iq,iquad,iqq;
422: PetscInt nli[num_z][2],indx[num_z];
423: PetscScalar dd,dl,zip,zipq,zz,bb,bbb,aij;
424: PetscScalar rquad[num_z][3],dlen[num_z],qdwt[3],add_term;
425: PetscErrorCode ierr;
427: /* initializing everything */
429: for(i=0; i < nz; i++)
430: {
431: nli[i][0] = 0;
432: nli[i][1] = 0;
433: indx[i] = 0;
434: rquad[i][0] = 0.0;
435: rquad[i][1] = 0.0;
436: rquad[i][2] = 0.0;
437: dlen[i] = 0.0;
438: }/*end for(i)*/
440: /* quadrature weights */
441: qdwt[0] = 1.0/6.0;
442: qdwt[1] = 4.0/6.0;
443: qdwt[2] = 1.0/6.0;
445: /* 1st and last nodes have Dirichlet boundary condition -
446: set indices there to -1 */
448: for(i=0; i < nz-1; i++)
449: {
450: indx[i]=i-1;
452: }/*end for(i)*/
453: indx[nz-1]=-1;
455: ipq = 0;
457: for(il=0; il < nz-1; il++)
458: {
459: ip = ipq;
460: ipq = ip+1;
461: zip = z[ip];
462: zipq = z[ipq];
463: dl = zipq-zip;
464: rquad[il][0] = zip;
465: rquad[il][1] = (0.5)*(zip+zipq);
466: rquad[il][2] = zipq;
467: dlen[il] = fabs(dl);
468: nli[il][0] = ip;
469: nli[il][1] = ipq;
471: }/*end for(il)*/
473: for(il=0; il < nz-1; il++){
474: for(iquad=0; iquad < 3; iquad++){
475: dd = (dlen[il])*(qdwt[iquad]);
476: zz = rquad[il][iquad];
478: for(iq=0; iq < 2; iq++){
479: ip = nli[il][iq];
480: bb = bspl(z,zz,il,iq,nli,1);
481: i = indx[ip];
482: if(i > -1){
483: for(iqq=0; iqq < 2; iqq++){
484: ipp = nli[il][iqq];
485: bbb = bspl(z,zz,il,iqq,nli,1);
486: j = indx[ipp];
487: aij = bb*bbb;
488: if(j > -1) {
489: add_term = aij*dd;
490: MatSetValue(obj->Amat,i,j,add_term,ADD_VALUES);
491: }/*endif*/
492: }/*end for(iqq)*/
493: }/*end if(i>0)*/
494: }/*end for(iq)*/
495: }/*end for(iquad)*/
496: }/*end for(il)*/
497: MatAssemblyBegin(obj->Amat,MAT_FINAL_ASSEMBLY);
498: MatAssemblyEnd(obj->Amat,MAT_FINAL_ASSEMBLY);
499: return;
500: }
502: /*---------------------------------------------------------
503: Function to fill the rhs vector with
504: By + g values ****
505: ---------------------------------------------------------*/
506: void rhs(AppCtx *obj,PetscScalar *y, PetscInt nz, PetscScalar *z, PetscReal t)
507: {
508: PetscInt i,j,js,je,jj;
509: PetscScalar val,g[num_z],btri[num_z][3],add_term;
510: PetscErrorCode ierr;
512: for(i=0; i < nz-2; i++){
513: for(j=0; j <= 2; j++){
514: btri[i][j]=0.0;
515: }
516: g[i] = 0.0;
517: }
519: /* call femBg to set the tri-diagonal b matrix and vector g */
520: femBg(btri,g,nz,z,t);
522: /* setting the entries of the right hand side vector */
523: for(i=0; i < nz-2; i++){
524: val = 0.0;
525: js = 0;
526: if(i == 0) js = 1;
527: je = 2;
528: if(i == nz-2) je = 1;
530: for(jj=js; jj <= je; jj++){
531: j = i+jj-1;
532: val += (btri[i][jj])*(y[j]);
533: }
534: add_term = val + g[i];
535: VecSetValue(obj->ksp_rhs,(PetscInt)i,(PetscScalar)add_term,INSERT_VALUES);
536: }
537: VecAssemblyBegin(obj->ksp_rhs);
538: VecAssemblyEnd(obj->ksp_rhs);
540: /* return to main driver function */
541: return;
542: }
544: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
545: %% Function to form the right hand side of the time-stepping problem. %%
546: %% -------------------------------------------------------------------------------------------%%
547: if (useAlhs):
548: globalout = By+g
549: else if (!useAlhs):
550: globalout = f(y,t)=Ainv(By+g),
551: in which the ksp solver to transform the problem A*ydot=By+g
552: to the problem ydot=f(y,t)=inv(A)*(By+g)
553: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
555: PetscErrorCode RHSfunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
556: {
558: AppCtx *obj = (AppCtx*)ctx;
559: PetscScalar *soln_ptr,soln[num_z-2];
560: PetscInt i,nz=obj->nz;
561: PetscReal time;
563: /* get the previous solution to compute updated system */
564: VecGetArray(globalin,&soln_ptr);
565: for(i=0;i < num_z-2;i++){
566: soln[i] = soln_ptr[i];
567: }
568: VecRestoreArray(globalin,&soln_ptr);
570: /* clear out the matrix and rhs for ksp to keep things straight */
571: VecSet(obj->ksp_rhs,(PetscScalar)0.0);
573: time = t;
574: /* get the updated system */
575: rhs(obj,soln,nz,obj->z,time); /* setup of the By+g rhs */
577: /* do a ksp solve to get the rhs for the ts problem */
578: if (obj->useAlhs){
579: /* ksp_sol = ksp_rhs */
580: VecCopy(obj->ksp_rhs,globalout);
581: } else {
582: /* ksp_sol = inv(Amat)*ksp_rhs */
583: Petsc_KSPSolve(obj);
584: VecCopy(obj->ksp_sol,globalout);
585: }
586: return 0;
587: }