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 10
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 PetscErrorCode femBg(PetscScalar[][3], PetscScalar *, PetscInt, PetscScalar *, PetscReal);
40: extern PetscErrorCode femA(AppCtx *, PetscInt, PetscScalar *);
41: extern PetscErrorCode rhs(AppCtx *, PetscScalar *, PetscInt, PetscScalar *, PetscReal);
42: extern PetscErrorCode RHSfunction(TS, PetscReal, Vec, Vec, void *);
44: int main(int argc, char **argv)
45: {
46: PetscInt i, m, nz, steps, max_steps, k, nphase = 1;
47: PetscScalar zInitial, zFinal, val, *z;
48: PetscReal stepsz[4], T, ftime;
49: TS ts;
50: SNES snes;
51: Mat Jmat;
52: AppCtx appctx; /* user-defined application context */
53: Vec init_sol; /* ts solution vector */
54: PetscMPIInt size;
57: PetscInitialize(&argc, &argv, (char *)0, help);
58: MPI_Comm_size(PETSC_COMM_WORLD, &size);
61: /* initializations */
62: zInitial = 0.0;
63: zFinal = 1.0;
64: nz = num_z;
65: m = nz - 2;
66: appctx.nz = nz;
67: max_steps = (PetscInt)10000;
69: appctx.m = m;
70: appctx.max_probsz = nz;
71: appctx.debug = PETSC_FALSE;
72: appctx.useAlhs = PETSC_FALSE;
74: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "", "");
75: PetscOptionsName("-debug", NULL, NULL, &appctx.debug);
76: PetscOptionsName("-useAlhs", NULL, NULL, &appctx.useAlhs);
77: PetscOptionsRangeInt("-nphase", NULL, NULL, nphase, &nphase, NULL, 1, 3);
78: PetscOptionsEnd();
79: T = 0.014 / nphase;
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, NULL, &appctx.Amat);
93: MatSetFromOptions(appctx.Amat);
94: MatSetUp(appctx.Amat);
95: /* set space grid points - interio points only! */
96: PetscMalloc1(nz + 1, &z);
97: for (i = 0; i < nz; i++) z[i] = (i) * ((zFinal - zInitial) / (nz - 1));
98: appctx.z = z;
99: femA(&appctx, nz, z);
101: /* create the jacobian matrix */
102: /*----------------------------*/
103: MatCreate(PETSC_COMM_WORLD, &Jmat);
104: MatSetSizes(Jmat, PETSC_DECIDE, PETSC_DECIDE, m, m);
105: MatSetFromOptions(Jmat);
106: MatSetUp(Jmat);
108: /* create working vectors for formulating rhs=inv(Alhs)*(Arhs*U + g) */
109: VecDuplicate(init_sol, &appctx.ksp_rhs);
110: VecDuplicate(init_sol, &appctx.ksp_sol);
112: /* set initial guess */
113: /*-------------------*/
114: for (i = 0; i < nz - 2; i++) {
115: val = exact(z[i + 1], 0.0);
116: VecSetValue(init_sol, i, (PetscScalar)val, INSERT_VALUES);
117: }
118: VecAssemblyBegin(init_sol);
119: VecAssemblyEnd(init_sol);
121: /*create a time-stepping context and set the problem type */
122: /*--------------------------------------------------------*/
123: TSCreate(PETSC_COMM_WORLD, &ts);
124: TSSetProblemType(ts, TS_NONLINEAR);
126: /* set time-step method */
127: TSSetType(ts, TSCN);
129: /* Set optional user-defined monitoring routine */
130: TSMonitorSet(ts, Monitor, &appctx, NULL);
131: /* set the right hand side of U_t = RHSfunction(U,t) */
132: TSSetRHSFunction(ts, NULL, (PetscErrorCode(*)(TS, PetscScalar, Vec, Vec, void *))RHSfunction, &appctx);
134: if (appctx.useAlhs) {
135: /* set the left hand side matrix of Amat*U_t = rhs(U,t) */
137: /* Note: this approach is incompatible with the finite differenced Jacobian set below because we can't restore the
138: * Alhs matrix without making a copy. Either finite difference the entire thing or use analytic Jacobians in both
139: * places.
140: */
141: TSSetIFunction(ts, NULL, TSComputeIFunctionLinear, &appctx);
142: TSSetIJacobian(ts, appctx.Amat, appctx.Amat, TSComputeIJacobianConstant, &appctx);
143: }
145: /* use petsc to compute the jacobian by finite differences */
146: TSGetSNES(ts, &snes);
147: SNESSetJacobian(snes, Jmat, Jmat, SNESComputeJacobianDefault, NULL);
149: /* get the command line options if there are any and set them */
150: TSSetFromOptions(ts);
152: #if defined(PETSC_HAVE_SUNDIALS2)
153: {
154: TSType type;
155: PetscBool sundialstype = PETSC_FALSE;
156: TSGetType(ts, &type);
157: PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &sundialstype);
159: }
160: #endif
161: /* Sets the initial solution */
162: TSSetSolution(ts, init_sol);
164: stepsz[0] = 1.0 / (2.0 * (nz - 1) * (nz - 1)); /* (mesh_size)^2/2.0 */
165: ftime = 0.0;
166: for (k = 0; k < nphase; k++) {
167: if (nphase > 1) PetscPrintf(PETSC_COMM_WORLD, "Phase %" PetscInt_FMT " initial time %g, stepsz %g, duration: %g\n", k, (double)ftime, (double)stepsz[k], (double)((k + 1) * T));
168: TSSetTime(ts, ftime);
169: TSSetTimeStep(ts, stepsz[k]);
170: TSSetMaxSteps(ts, max_steps);
171: TSSetMaxTime(ts, (k + 1) * T);
172: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
174: /* loop over time steps */
175: /*----------------------*/
176: TSSolve(ts, init_sol);
177: TSGetSolveTime(ts, &ftime);
178: TSGetStepNumber(ts, &steps);
179: stepsz[k + 1] = stepsz[k] * 1.5; /* change step size for the next phase */
180: }
182: /* free space */
183: TSDestroy(&ts);
184: MatDestroy(&appctx.Amat);
185: MatDestroy(&Jmat);
186: VecDestroy(&appctx.ksp_rhs);
187: VecDestroy(&appctx.ksp_sol);
188: VecDestroy(&init_sol);
189: VecDestroy(&appctx.solution);
190: PetscFree(z);
192: PetscFinalize();
193: return 0;
194: }
196: /*------------------------------------------------------------------------
197: Set exact solution
198: u(z,t) = sin(6*PI*z)*exp(-36.*PI*PI*t) + 3.*sin(2*PI*z)*exp(-4.*PI*PI*t)
199: --------------------------------------------------------------------------*/
200: PetscScalar exact(PetscScalar z, PetscReal t)
201: {
202: PetscScalar val, ex1, ex2;
204: ex1 = PetscExpReal(-36. * PETSC_PI * PETSC_PI * t);
205: ex2 = PetscExpReal(-4. * PETSC_PI * PETSC_PI * t);
206: val = PetscSinScalar(6 * PETSC_PI * z) * ex1 + 3. * PetscSinScalar(2 * PETSC_PI * z) * ex2;
207: return val;
208: }
210: /*
211: Monitor - User-provided routine to monitor the solution computed at
212: each timestep. This example plots the solution and computes the
213: error in two different norms.
215: Input Parameters:
216: ts - the timestep context
217: step - the count of the current step (with 0 meaning the
218: initial condition)
219: time - the current time
220: u - the solution at this timestep
221: ctx - the user-provided context for this monitoring routine.
222: In this case we use the application context which contains
223: information about the problem size, workspace and the exact
224: solution.
225: */
226: PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec u, void *ctx)
227: {
228: AppCtx *appctx = (AppCtx *)ctx;
229: PetscInt i, m = appctx->m;
230: PetscReal norm_2, norm_max, h = 1.0 / (m + 1);
231: PetscScalar *u_exact;
233: /* Compute the exact solution */
234: VecGetArrayWrite(appctx->solution, &u_exact);
235: for (i = 0; i < m; i++) u_exact[i] = exact(appctx->z[i + 1], time);
236: VecRestoreArrayWrite(appctx->solution, &u_exact);
238: /* Print debugging information if desired */
239: if (appctx->debug) {
240: PetscPrintf(PETSC_COMM_SELF, "Computed solution vector at time %g\n", (double)time);
241: VecView(u, PETSC_VIEWER_STDOUT_SELF);
242: PetscPrintf(PETSC_COMM_SELF, "Exact solution vector\n");
243: VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF);
244: }
246: /* Compute the 2-norm and max-norm of the error */
247: VecAXPY(appctx->solution, -1.0, u);
248: VecNorm(appctx->solution, NORM_2, &norm_2);
250: norm_2 = PetscSqrtReal(h) * norm_2;
251: VecNorm(appctx->solution, NORM_MAX, &norm_max);
252: PetscPrintf(PETSC_COMM_SELF, "Timestep %" PetscInt_FMT ": time = %g, 2-norm error = %6.4f, max norm error = %6.4f\n", step, (double)time, (double)norm_2, (double)norm_max);
254: /*
255: Print debugging information if desired
256: */
257: if (appctx->debug) {
258: PetscPrintf(PETSC_COMM_SELF, "Error vector\n");
259: VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF);
260: }
261: return 0;
262: }
264: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
265: Function to solve a linear system using KSP
266: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
268: PetscErrorCode Petsc_KSPSolve(AppCtx *obj)
269: {
270: KSP ksp;
271: PC pc;
273: /*create the ksp context and set the operators,that is, associate the system matrix with it*/
274: KSPCreate(PETSC_COMM_WORLD, &ksp);
275: KSPSetOperators(ksp, obj->Amat, obj->Amat);
277: /*get the preconditioner context, set its type and the tolerances*/
278: KSPGetPC(ksp, &pc);
279: PCSetType(pc, PCLU);
280: KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
282: /*get the command line options if there are any and set them*/
283: KSPSetFromOptions(ksp);
285: /*get the linear system (ksp) solve*/
286: KSPSolve(ksp, obj->ksp_rhs, obj->ksp_sol);
288: KSPDestroy(&ksp);
289: return 0;
290: }
292: /***********************************************************************
293: Function to return value of basis function or derivative of basis function.
294: ***********************************************************************
296: Arguments:
297: x = array of xpoints or nodal values
298: xx = point at which the basis function is to be
299: evaluated.
300: il = interval containing xx.
301: iq = indicates which of the two basis functions in
302: interval intrvl should be used
303: nll = array containing the endpoints of each interval.
304: id = If id ~= 2, the value of the basis function
305: is calculated; if id = 2, the value of the
306: derivative of the basis function is returned.
307: ***********************************************************************/
309: PetscScalar bspl(PetscScalar *x, PetscScalar xx, PetscInt il, PetscInt iq, PetscInt nll[][2], PetscInt id)
310: {
311: PetscScalar x1, x2, bfcn;
312: PetscInt i1, i2, iq1, iq2;
314: /* Determine which basis function in interval intrvl is to be used in */
315: iq1 = iq;
316: if (iq1 == 0) iq2 = 1;
317: else iq2 = 0;
319: /* Determine endpoint of the interval intrvl */
320: i1 = nll[il][iq1];
321: i2 = nll[il][iq2];
323: /* Determine nodal values at the endpoints of the interval intrvl */
324: x1 = x[i1];
325: x2 = x[i2];
327: /* Evaluate basis function */
328: if (id == 2) bfcn = (1.0) / (x1 - x2);
329: else bfcn = (xx - x2) / (x1 - x2);
330: return bfcn;
331: }
333: /*---------------------------------------------------------
334: Function called by rhs function to get B and g
335: ---------------------------------------------------------*/
336: PetscErrorCode 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, b_z, 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++) indx[i] = i - 1;
363: indx[nz - 1] = -1;
365: ipq = 0;
366: for (il = 0; il < nz - 1; il++) {
367: ip = ipq;
368: ipq = ip + 1;
369: zip = z[ip];
370: zipq = z[ipq];
371: dl = zipq - zip;
372: zquad[il][0] = zip;
373: zquad[il][1] = (0.5) * (zip + zipq);
374: zquad[il][2] = zipq;
375: dlen[il] = PetscAbsScalar(dl);
376: nli[il][0] = ip;
377: nli[il][1] = ipq;
378: }
380: for (il = 0; il < nz - 1; il++) {
381: for (iquad = 0; iquad < 3; iquad++) {
382: dd = (dlen[il]) * (qdwt[iquad]);
383: zz = zquad[il][iquad];
385: for (iq = 0; iq < 2; iq++) {
386: ip = nli[il][iq];
387: b_z = bspl(z, zz, il, iq, nli, 2);
388: i = indx[ip];
390: if (i > -1) {
391: for (iqq = 0; iqq < 2; iqq++) {
392: ipp = nli[il][iqq];
393: bb_z = bspl(z, zz, il, iqq, nli, 2);
394: j = indx[ipp];
395: bij = -b_z * bb_z;
397: if (j > -1) {
398: jj = 1 + j - i;
399: btri[i][jj] += bij * dd;
400: } else {
401: f[i] += bij * dd * exact(z[ipp], t);
402: /* f[i] += 0.0; */
403: /* if (il==0 && j==-1) { */
404: /* f[i] += bij*dd*exact(zz,t); */
405: /* }*/ /*end if*/
406: } /*end else*/
407: } /*end for (iqq)*/
408: } /*end if (i>0)*/
409: } /*end for (iq)*/
410: } /*end for (iquad)*/
411: } /*end for (il)*/
412: return 0;
413: }
415: PetscErrorCode femA(AppCtx *obj, PetscInt nz, PetscScalar *z)
416: {
417: PetscInt i, j, il, ip, ipp, ipq, iq, iquad, iqq;
418: PetscInt nli[num_z][2], indx[num_z];
419: PetscScalar dd, dl, zip, zipq, zz, bb, bbb, aij;
420: PetscScalar rquad[num_z][3], dlen[num_z], qdwt[3], add_term;
422: /* initializing everything */
423: for (i = 0; i < nz; i++) {
424: nli[i][0] = 0;
425: nli[i][1] = 0;
426: indx[i] = 0;
427: rquad[i][0] = 0.0;
428: rquad[i][1] = 0.0;
429: rquad[i][2] = 0.0;
430: dlen[i] = 0.0;
431: } /*end for (i)*/
433: /* quadrature weights */
434: qdwt[0] = 1.0 / 6.0;
435: qdwt[1] = 4.0 / 6.0;
436: qdwt[2] = 1.0 / 6.0;
438: /* 1st and last nodes have Dirichlet boundary condition -
439: set indices there to -1 */
441: for (i = 0; i < nz - 1; i++) indx[i] = i - 1;
442: indx[nz - 1] = -1;
444: ipq = 0;
446: for (il = 0; il < nz - 1; il++) {
447: ip = ipq;
448: ipq = ip + 1;
449: zip = z[ip];
450: zipq = z[ipq];
451: dl = zipq - zip;
452: rquad[il][0] = zip;
453: rquad[il][1] = (0.5) * (zip + zipq);
454: rquad[il][2] = zipq;
455: dlen[il] = PetscAbsScalar(dl);
456: nli[il][0] = ip;
457: nli[il][1] = ipq;
458: } /*end for (il)*/
460: for (il = 0; il < nz - 1; il++) {
461: for (iquad = 0; iquad < 3; iquad++) {
462: dd = (dlen[il]) * (qdwt[iquad]);
463: zz = rquad[il][iquad];
465: for (iq = 0; iq < 2; iq++) {
466: ip = nli[il][iq];
467: bb = bspl(z, zz, il, iq, nli, 1);
468: i = indx[ip];
469: if (i > -1) {
470: for (iqq = 0; iqq < 2; iqq++) {
471: ipp = nli[il][iqq];
472: bbb = bspl(z, zz, il, iqq, nli, 1);
473: j = indx[ipp];
474: aij = bb * bbb;
475: if (j > -1) {
476: add_term = aij * dd;
477: MatSetValue(obj->Amat, i, j, add_term, ADD_VALUES);
478: } /*endif*/
479: } /*end for (iqq)*/
480: } /*end if (i>0)*/
481: } /*end for (iq)*/
482: } /*end for (iquad)*/
483: } /*end for (il)*/
484: MatAssemblyBegin(obj->Amat, MAT_FINAL_ASSEMBLY);
485: MatAssemblyEnd(obj->Amat, MAT_FINAL_ASSEMBLY);
486: return 0;
487: }
489: /*---------------------------------------------------------
490: Function to fill the rhs vector with
491: By + g values ****
492: ---------------------------------------------------------*/
493: PetscErrorCode rhs(AppCtx *obj, PetscScalar *y, PetscInt nz, PetscScalar *z, PetscReal t)
494: {
495: PetscInt i, j, js, je, jj;
496: PetscScalar val, g[num_z], btri[num_z][3], add_term;
498: for (i = 0; i < nz - 2; i++) {
499: for (j = 0; j <= 2; j++) btri[i][j] = 0.0;
500: g[i] = 0.0;
501: }
503: /* call femBg to set the tri-diagonal b matrix and vector g */
504: femBg(btri, g, nz, z, t);
506: /* setting the entries of the right hand side vector */
507: for (i = 0; i < nz - 2; i++) {
508: val = 0.0;
509: js = 0;
510: if (i == 0) js = 1;
511: je = 2;
512: if (i == nz - 2) je = 1;
514: for (jj = js; jj <= je; jj++) {
515: j = i + jj - 1;
516: val += (btri[i][jj]) * (y[j]);
517: }
518: add_term = val + g[i];
519: VecSetValue(obj->ksp_rhs, (PetscInt)i, (PetscScalar)add_term, INSERT_VALUES);
520: }
521: VecAssemblyBegin(obj->ksp_rhs);
522: VecAssemblyEnd(obj->ksp_rhs);
523: return 0;
524: }
526: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
527: %% Function to form the right hand side of the time-stepping problem. %%
528: %% -------------------------------------------------------------------------------------------%%
529: if (useAlhs):
530: globalout = By+g
531: else if (!useAlhs):
532: globalout = f(y,t)=Ainv(By+g),
533: in which the ksp solver to transform the problem A*ydot=By+g
534: to the problem ydot=f(y,t)=inv(A)*(By+g)
535: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
537: PetscErrorCode RHSfunction(TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx)
538: {
539: AppCtx *obj = (AppCtx *)ctx;
540: PetscScalar soln[num_z];
541: const PetscScalar *soln_ptr;
542: PetscInt i, nz = obj->nz;
543: PetscReal time;
545: /* get the previous solution to compute updated system */
546: VecGetArrayRead(globalin, &soln_ptr);
547: for (i = 0; i < num_z - 2; i++) soln[i] = soln_ptr[i];
548: VecRestoreArrayRead(globalin, &soln_ptr);
549: soln[num_z - 1] = 0.0;
550: soln[num_z - 2] = 0.0;
552: /* clear out the matrix and rhs for ksp to keep things straight */
553: VecSet(obj->ksp_rhs, (PetscScalar)0.0);
555: time = t;
556: /* get the updated system */
557: rhs(obj, soln, nz, obj->z, time); /* setup of the By+g rhs */
559: /* do a ksp solve to get the rhs for the ts problem */
560: if (obj->useAlhs) {
561: /* ksp_sol = ksp_rhs */
562: VecCopy(obj->ksp_rhs, globalout);
563: } else {
564: /* ksp_sol = inv(Amat)*ksp_rhs */
565: Petsc_KSPSolve(obj);
566: VecCopy(obj->ksp_sol, globalout);
567: }
568: return 0;
569: }
571: /*TEST
573: build:
574: requires: !complex
576: test:
577: suffix: euler
578: output_file: output/ex3.out
580: test:
581: suffix: 2
582: args: -useAlhs
583: output_file: output/ex3.out
584: TODO: Broken because SNESComputeJacobianDefault is incompatible with TSComputeIJacobianConstant
586: TEST*/