Actual source code: ex3.c
1: static char help[] = "Solves 1D heat equation with FEM formulation.\n\
2: Input arguments are\n\
3: -useAlhs: solve Alhs*U' = (Arhs*U + g) \n\
4: otherwise, solve U' = inv(Alhs)*(Arhs*U + g) \n\n";
6: /*--------------------------------------------------------------------------
7: Solves 1D heat equation U_t = U_xx with FEM formulation:
8: Alhs*U' = rhs (= Arhs*U + g)
9: We thank Chris Cox <clcox@clemson.edu> for contributing the original code
10: ----------------------------------------------------------------------------*/
12: #include <petscksp.h>
13: #include <petscts.h>
15: /* special variable - max size of all arrays */
16: #define num_z 10
18: /*
19: User-defined application context - contains data needed by the
20: application-provided call-back routines.
21: */
22: typedef struct {
23: Mat Amat; /* left hand side matrix */
24: Vec ksp_rhs, ksp_sol; /* working vectors for formulating inv(Alhs)*(Arhs*U+g) */
25: PetscInt max_probsz; /* max size of the problem */
26: PetscBool useAlhs; /* flag (1 indicates solving Alhs*U' = Arhs*U+g */
27: PetscInt nz; /* total number of grid points */
28: PetscInt m; /* total number of interio grid points */
29: Vec solution; /* global exact ts solution vector */
30: PetscScalar *z; /* array of grid points */
31: PetscBool debug; /* flag (1 indicates activation of debugging printouts) */
32: } AppCtx;
34: extern PetscScalar exact(PetscScalar, PetscReal);
35: extern PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *);
36: extern PetscErrorCode Petsc_KSPSolve(AppCtx *);
37: extern PetscScalar bspl(PetscScalar *, PetscScalar, PetscInt, PetscInt, PetscInt[][2], PetscInt);
38: extern PetscErrorCode femBg(PetscScalar[][3], PetscScalar *, PetscInt, PetscScalar *, PetscReal);
39: extern PetscErrorCode femA(AppCtx *, PetscInt, PetscScalar *);
40: extern PetscErrorCode rhs(AppCtx *, PetscScalar *, PetscInt, PetscScalar *, PetscReal);
41: extern PetscErrorCode RHSfunction(TS, PetscReal, Vec, Vec, void *);
43: int main(int argc, char **argv)
44: {
45: PetscInt i, m, nz, steps, max_steps, k, nphase = 1;
46: PetscScalar zInitial, zFinal, val, *z;
47: PetscReal stepsz[4], T, ftime;
48: TS ts;
49: SNES snes;
50: Mat Jmat;
51: AppCtx appctx; /* user-defined application context */
52: Vec init_sol; /* ts solution vector */
53: PetscMPIInt size;
55: PetscFunctionBeginUser;
56: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
57: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
58: PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only");
60: /* initializations */
61: zInitial = 0.0;
62: zFinal = 1.0;
63: nz = num_z;
64: m = nz - 2;
65: appctx.nz = nz;
66: max_steps = (PetscInt)10000;
68: appctx.m = m;
69: appctx.max_probsz = nz;
70: appctx.debug = PETSC_FALSE;
71: appctx.useAlhs = PETSC_FALSE;
73: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "", "");
74: PetscCall(PetscOptionsName("-debug", NULL, NULL, &appctx.debug));
75: PetscCall(PetscOptionsName("-useAlhs", NULL, NULL, &appctx.useAlhs));
76: PetscCall(PetscOptionsRangeInt("-nphase", NULL, NULL, nphase, &nphase, NULL, 1, 3));
77: PetscOptionsEnd();
78: T = 0.014 / nphase;
80: /* create vector to hold ts solution */
81: /*-----------------------------------*/
82: PetscCall(VecCreate(PETSC_COMM_WORLD, &init_sol));
83: PetscCall(VecSetSizes(init_sol, PETSC_DECIDE, m));
84: PetscCall(VecSetFromOptions(init_sol));
86: /* create vector to hold true ts soln for comparison */
87: PetscCall(VecDuplicate(init_sol, &appctx.solution));
89: /* create LHS matrix Amat */
90: /*------------------------*/
91: PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, m, m, 3, NULL, &appctx.Amat));
92: PetscCall(MatSetFromOptions(appctx.Amat));
93: PetscCall(MatSetUp(appctx.Amat));
94: /* set space grid points - interio points only! */
95: PetscCall(PetscMalloc1(nz + 1, &z));
96: for (i = 0; i < nz; i++) z[i] = (i) * ((zFinal - zInitial) / (nz - 1));
97: appctx.z = z;
98: PetscCall(femA(&appctx, nz, z));
100: /* create the jacobian matrix */
101: /*----------------------------*/
102: PetscCall(MatCreate(PETSC_COMM_WORLD, &Jmat));
103: PetscCall(MatSetSizes(Jmat, PETSC_DECIDE, PETSC_DECIDE, m, m));
104: PetscCall(MatSetFromOptions(Jmat));
105: PetscCall(MatSetUp(Jmat));
107: /* create working vectors for formulating rhs=inv(Alhs)*(Arhs*U + g) */
108: PetscCall(VecDuplicate(init_sol, &appctx.ksp_rhs));
109: PetscCall(VecDuplicate(init_sol, &appctx.ksp_sol));
111: /* set initial guess */
112: /*-------------------*/
113: for (i = 0; i < nz - 2; i++) {
114: val = exact(z[i + 1], 0.0);
115: PetscCall(VecSetValue(init_sol, i, (PetscScalar)val, INSERT_VALUES));
116: }
117: PetscCall(VecAssemblyBegin(init_sol));
118: PetscCall(VecAssemblyEnd(init_sol));
120: /*create a time-stepping context and set the problem type */
121: /*--------------------------------------------------------*/
122: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
123: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
125: /* set time-step method */
126: PetscCall(TSSetType(ts, TSCN));
128: /* Set optional user-defined monitoring routine */
129: PetscCall(TSMonitorSet(ts, Monitor, &appctx, NULL));
130: /* set the right-hand side of U_t = RHSfunction(U,t) */
131: PetscCall(TSSetRHSFunction(ts, NULL, (PetscErrorCode (*)(TS, PetscScalar, Vec, Vec, void *))RHSfunction, &appctx));
133: if (appctx.useAlhs) {
134: /* set the left hand side matrix of Amat*U_t = rhs(U,t) */
136: /* Note: this approach is incompatible with the finite differenced Jacobian set below because we can't restore the
137: * Alhs matrix without making a copy. Either finite difference the entire thing or use analytic Jacobians in both
138: * places.
139: */
140: PetscCall(TSSetIFunction(ts, NULL, TSComputeIFunctionLinear, &appctx));
141: PetscCall(TSSetIJacobian(ts, appctx.Amat, appctx.Amat, TSComputeIJacobianConstant, &appctx));
142: }
144: /* use petsc to compute the jacobian by finite differences */
145: PetscCall(TSGetSNES(ts, &snes));
146: PetscCall(SNESSetJacobian(snes, Jmat, Jmat, SNESComputeJacobianDefault, NULL));
148: /* get the command line options if there are any and set them */
149: PetscCall(TSSetFromOptions(ts));
151: #if defined(PETSC_HAVE_SUNDIALS2)
152: {
153: TSType type;
154: PetscBool sundialstype = PETSC_FALSE;
155: PetscCall(TSGetType(ts, &type));
156: PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &sundialstype));
157: PetscCheck(!sundialstype || !appctx.useAlhs, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot use Alhs formulation for TSSUNDIALS type");
158: }
159: #endif
160: /* Sets the initial solution */
161: PetscCall(TSSetSolution(ts, init_sol));
163: stepsz[0] = 1.0 / (2.0 * (nz - 1) * (nz - 1)); /* (mesh_size)^2/2.0 */
164: ftime = 0.0;
165: for (k = 0; k < nphase; k++) {
166: if (nphase > 1) PetscCall(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)));
167: PetscCall(TSSetTime(ts, ftime));
168: PetscCall(TSSetTimeStep(ts, stepsz[k]));
169: PetscCall(TSSetMaxSteps(ts, max_steps));
170: PetscCall(TSSetMaxTime(ts, (k + 1) * T));
171: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
173: /* loop over time steps */
174: /*----------------------*/
175: PetscCall(TSSolve(ts, init_sol));
176: PetscCall(TSGetSolveTime(ts, &ftime));
177: PetscCall(TSGetStepNumber(ts, &steps));
178: stepsz[k + 1] = stepsz[k] * 1.5; /* change step size for the next phase */
179: }
181: /* free space */
182: PetscCall(TSDestroy(&ts));
183: PetscCall(MatDestroy(&appctx.Amat));
184: PetscCall(MatDestroy(&Jmat));
185: PetscCall(VecDestroy(&appctx.ksp_rhs));
186: PetscCall(VecDestroy(&appctx.ksp_sol));
187: PetscCall(VecDestroy(&init_sol));
188: PetscCall(VecDestroy(&appctx.solution));
189: PetscCall(PetscFree(z));
191: PetscCall(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 = PetscExpReal(-36. * PETSC_PI * PETSC_PI * t);
204: ex2 = PetscExpReal(-4. * PETSC_PI * PETSC_PI * t);
205: val = PetscSinScalar(6 * PETSC_PI * z) * ex1 + 3. * PetscSinScalar(2 * PETSC_PI * z) * ex2;
206: return val;
207: }
209: /*
210: Monitor - User-provided routine to monitor the solution computed at
211: each timestep. This example plots the solution and computes the
212: error in two different norms.
214: Input Parameters:
215: ts - the timestep context
216: step - the count of the current step (with 0 meaning the
217: initial condition)
218: time - the current time
219: u - the solution at this timestep
220: ctx - the user-provided context for this monitoring routine.
221: In this case we use the application context which contains
222: information about the problem size, workspace and the exact
223: solution.
224: */
225: PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec u, void *ctx)
226: {
227: AppCtx *appctx = (AppCtx *)ctx;
228: PetscInt i, m = appctx->m;
229: PetscReal norm_2, norm_max, h = 1.0 / (m + 1);
230: PetscScalar *u_exact;
232: PetscFunctionBeginUser;
233: /* Compute the exact solution */
234: PetscCall(VecGetArrayWrite(appctx->solution, &u_exact));
235: for (i = 0; i < m; i++) u_exact[i] = exact(appctx->z[i + 1], time);
236: PetscCall(VecRestoreArrayWrite(appctx->solution, &u_exact));
238: /* Print debugging information if desired */
239: if (appctx->debug) {
240: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Computed solution vector at time %g\n", (double)time));
241: PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF));
242: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Exact solution vector\n"));
243: PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
244: }
246: /* Compute the 2-norm and max-norm of the error */
247: PetscCall(VecAXPY(appctx->solution, -1.0, u));
248: PetscCall(VecNorm(appctx->solution, NORM_2, &norm_2));
250: norm_2 = PetscSqrtReal(h) * norm_2;
251: PetscCall(VecNorm(appctx->solution, NORM_MAX, &norm_max));
252: PetscCall(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: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error vector\n"));
259: PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
260: }
261: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBeginUser;
274: /*create the ksp context and set the operators,that is, associate the system matrix with it*/
275: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
276: PetscCall(KSPSetOperators(ksp, obj->Amat, obj->Amat));
278: /*get the preconditioner context, set its type and the tolerances*/
279: PetscCall(KSPGetPC(ksp, &pc));
280: PetscCall(PCSetType(pc, PCLU));
281: PetscCall(KSPSetTolerances(ksp, 1.e-7, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
283: /*get the command line options if there are any and set them*/
284: PetscCall(KSPSetFromOptions(ksp));
286: /*get the linear system (ksp) solve*/
287: PetscCall(KSPSolve(ksp, obj->ksp_rhs, obj->ksp_sol));
289: PetscCall(KSPDestroy(&ksp));
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: /***********************************************************************
294: Function to return value of basis function or derivative of basis function.
295: ***********************************************************************
297: Arguments:
298: x = array of xpoints or nodal values
299: xx = point at which the basis function is to be
300: evaluated.
301: il = interval containing xx.
302: iq = indicates which of the two basis functions in
303: interval intrvl should be used
304: nll = array containing the endpoints of each interval.
305: id = If id ~= 2, the value of the basis function
306: is calculated; if id = 2, the value of the
307: derivative of the basis function is returned.
308: ***********************************************************************/
310: PetscScalar bspl(PetscScalar *x, PetscScalar xx, PetscInt il, PetscInt iq, PetscInt nll[][2], PetscInt id)
311: {
312: PetscScalar x1, x2, bfcn;
313: PetscInt i1, i2, iq1, iq2;
315: /* Determine which basis function in interval intrvl is to be used in */
316: iq1 = iq;
317: if (iq1 == 0) iq2 = 1;
318: else iq2 = 0;
320: /* Determine endpoint of the interval intrvl */
321: i1 = nll[il][iq1];
322: i2 = nll[il][iq2];
324: /* Determine nodal values at the endpoints of the interval intrvl */
325: x1 = x[i1];
326: x2 = x[i2];
328: /* Evaluate basis function */
329: if (id == 2) bfcn = (1.0) / (x1 - x2);
330: else bfcn = (xx - x2) / (x1 - x2);
331: return bfcn;
332: }
334: /*---------------------------------------------------------
335: Function called by rhs function to get B and g
336: ---------------------------------------------------------*/
337: PetscErrorCode femBg(PetscScalar btri[][3], PetscScalar *f, PetscInt nz, PetscScalar *z, PetscReal t)
338: {
339: PetscInt i, j, jj, il, ip, ipp, ipq, iq, iquad, iqq;
340: PetscInt nli[num_z][2], indx[num_z];
341: PetscScalar dd, dl, zip, zipq, zz, b_z, bb_z, bij;
342: PetscScalar zquad[num_z][3], dlen[num_z], qdwt[3];
344: PetscFunctionBeginUser;
345: /* initializing everything - btri and f are initialized in rhs.c */
346: for (i = 0; i < nz; i++) {
347: nli[i][0] = 0;
348: nli[i][1] = 0;
349: indx[i] = 0;
350: zquad[i][0] = 0.0;
351: zquad[i][1] = 0.0;
352: zquad[i][2] = 0.0;
353: dlen[i] = 0.0;
354: } /*end for (i)*/
356: /* quadrature weights */
357: qdwt[0] = 1.0 / 6.0;
358: qdwt[1] = 4.0 / 6.0;
359: qdwt[2] = 1.0 / 6.0;
361: /* 1st and last nodes have Dirichlet boundary condition -
362: set indices there to -1 */
364: for (i = 0; i < nz - 1; i++) indx[i] = i - 1;
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] = PetscAbsScalar(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: b_z = bspl(z, zz, il, iq, nli, 2);
390: i = indx[ip];
392: if (i > -1) {
393: for (iqq = 0; iqq < 2; iqq++) {
394: ipp = nli[il][iqq];
395: bb_z = bspl(z, zz, il, iqq, nli, 2);
396: j = indx[ipp];
397: bij = -b_z * bb_z;
399: if (j > -1) {
400: jj = 1 + j - i;
401: btri[i][jj] += bij * dd;
402: } else {
403: f[i] += bij * dd * exact(z[ipp], t);
404: /* f[i] += 0.0; */
405: /* if (il==0 && j==-1) { */
406: /* f[i] += bij*dd*exact(zz,t); */
407: /* }*/ /*end if*/
408: } /*end else*/
409: } /*end for (iqq)*/
410: } /*end if (i>0)*/
411: } /*end for (iq)*/
412: } /*end for (iquad)*/
413: } /*end for (il)*/
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: PetscErrorCode femA(AppCtx *obj, PetscInt nz, PetscScalar *z)
418: {
419: PetscInt i, j, il, ip, ipp, ipq, iq, iquad, iqq;
420: PetscInt nli[num_z][2], indx[num_z];
421: PetscScalar dd, dl, zip, zipq, zz, bb, bbb, aij;
422: PetscScalar rquad[num_z][3], dlen[num_z], qdwt[3], add_term;
424: PetscFunctionBeginUser;
425: /* initializing everything */
426: for (i = 0; i < nz; i++) {
427: nli[i][0] = 0;
428: nli[i][1] = 0;
429: indx[i] = 0;
430: rquad[i][0] = 0.0;
431: rquad[i][1] = 0.0;
432: rquad[i][2] = 0.0;
433: dlen[i] = 0.0;
434: } /*end for (i)*/
436: /* quadrature weights */
437: qdwt[0] = 1.0 / 6.0;
438: qdwt[1] = 4.0 / 6.0;
439: qdwt[2] = 1.0 / 6.0;
441: /* 1st and last nodes have Dirichlet boundary condition -
442: set indices there to -1 */
444: for (i = 0; i < nz - 1; i++) indx[i] = i - 1;
445: indx[nz - 1] = -1;
447: ipq = 0;
449: for (il = 0; il < nz - 1; il++) {
450: ip = ipq;
451: ipq = ip + 1;
452: zip = z[ip];
453: zipq = z[ipq];
454: dl = zipq - zip;
455: rquad[il][0] = zip;
456: rquad[il][1] = (0.5) * (zip + zipq);
457: rquad[il][2] = zipq;
458: dlen[il] = PetscAbsScalar(dl);
459: nli[il][0] = ip;
460: nli[il][1] = ipq;
461: } /*end for (il)*/
463: for (il = 0; il < nz - 1; il++) {
464: for (iquad = 0; iquad < 3; iquad++) {
465: dd = (dlen[il]) * (qdwt[iquad]);
466: zz = rquad[il][iquad];
468: for (iq = 0; iq < 2; iq++) {
469: ip = nli[il][iq];
470: bb = bspl(z, zz, il, iq, nli, 1);
471: i = indx[ip];
472: if (i > -1) {
473: for (iqq = 0; iqq < 2; iqq++) {
474: ipp = nli[il][iqq];
475: bbb = bspl(z, zz, il, iqq, nli, 1);
476: j = indx[ipp];
477: aij = bb * bbb;
478: if (j > -1) {
479: add_term = aij * dd;
480: PetscCall(MatSetValue(obj->Amat, i, j, add_term, ADD_VALUES));
481: } /*endif*/
482: } /*end for (iqq)*/
483: } /*end if (i>0)*/
484: } /*end for (iq)*/
485: } /*end for (iquad)*/
486: } /*end for (il)*/
487: PetscCall(MatAssemblyBegin(obj->Amat, MAT_FINAL_ASSEMBLY));
488: PetscCall(MatAssemblyEnd(obj->Amat, MAT_FINAL_ASSEMBLY));
489: PetscFunctionReturn(PETSC_SUCCESS);
490: }
492: /*---------------------------------------------------------
493: Function to fill the rhs vector with
494: By + g values ****
495: ---------------------------------------------------------*/
496: PetscErrorCode rhs(AppCtx *obj, PetscScalar *y, PetscInt nz, PetscScalar *z, PetscReal t)
497: {
498: PetscInt i, j, js, je, jj;
499: PetscScalar val, g[num_z], btri[num_z][3], add_term;
501: PetscFunctionBeginUser;
502: for (i = 0; i < nz - 2; i++) {
503: for (j = 0; j <= 2; j++) btri[i][j] = 0.0;
504: g[i] = 0.0;
505: }
507: /* call femBg to set the tri-diagonal b matrix and vector g */
508: PetscCall(femBg(btri, g, nz, z, t));
510: /* setting the entries of the right-hand side vector */
511: for (i = 0; i < nz - 2; i++) {
512: val = 0.0;
513: js = 0;
514: if (i == 0) js = 1;
515: je = 2;
516: if (i == nz - 2) je = 1;
518: for (jj = js; jj <= je; jj++) {
519: j = i + jj - 1;
520: val += (btri[i][jj]) * (y[j]);
521: }
522: add_term = val + g[i];
523: PetscCall(VecSetValue(obj->ksp_rhs, (PetscInt)i, (PetscScalar)add_term, INSERT_VALUES));
524: }
525: PetscCall(VecAssemblyBegin(obj->ksp_rhs));
526: PetscCall(VecAssemblyEnd(obj->ksp_rhs));
527: PetscFunctionReturn(PETSC_SUCCESS);
528: }
530: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
531: %% Function to form the right-hand side of the time-stepping problem. %%
532: %% -------------------------------------------------------------------------------------------%%
533: if (useAlhs):
534: globalout = By+g
535: else if (!useAlhs):
536: globalout = f(y,t)=Ainv(By+g),
537: in which the ksp solver to transform the problem A*ydot=By+g
538: to the problem ydot=f(y,t)=inv(A)*(By+g)
539: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
541: PetscErrorCode RHSfunction(TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx)
542: {
543: AppCtx *obj = (AppCtx *)ctx;
544: PetscScalar soln[num_z];
545: const PetscScalar *soln_ptr;
546: PetscInt i, nz = obj->nz;
547: PetscReal time;
549: PetscFunctionBeginUser;
550: /* get the previous solution to compute updated system */
551: PetscCall(VecGetArrayRead(globalin, &soln_ptr));
552: for (i = 0; i < num_z - 2; i++) soln[i] = soln_ptr[i];
553: PetscCall(VecRestoreArrayRead(globalin, &soln_ptr));
554: soln[num_z - 1] = 0.0;
555: soln[num_z - 2] = 0.0;
557: /* clear out the matrix and rhs for ksp to keep things straight */
558: PetscCall(VecSet(obj->ksp_rhs, (PetscScalar)0.0));
560: time = t;
561: /* get the updated system */
562: PetscCall(rhs(obj, soln, nz, obj->z, time)); /* setup of the By+g rhs */
564: /* do a ksp solve to get the rhs for the ts problem */
565: if (obj->useAlhs) {
566: /* ksp_sol = ksp_rhs */
567: PetscCall(VecCopy(obj->ksp_rhs, globalout));
568: } else {
569: /* ksp_sol = inv(Amat)*ksp_rhs */
570: PetscCall(Petsc_KSPSolve(obj));
571: PetscCall(VecCopy(obj->ksp_sol, globalout));
572: }
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: /*TEST
578: build:
579: requires: !complex
581: test:
582: suffix: euler
583: output_file: output/ex3.out
585: test:
586: suffix: 2
587: args: -useAlhs
588: output_file: output/ex3.out
589: TODO: Broken because SNESComputeJacobianDefault is incompatible with TSComputeIJacobianConstant
591: TEST*/