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;

 56:   PetscFunctionBeginUser;
 57:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 58:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 59:   PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only");

 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:   PetscCall(PetscOptionsName("-debug", NULL, NULL, &appctx.debug));
 76:   PetscCall(PetscOptionsName("-useAlhs", NULL, NULL, &appctx.useAlhs));
 77:   PetscCall(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:   PetscCall(VecCreate(PETSC_COMM_WORLD, &init_sol));
 84:   PetscCall(VecSetSizes(init_sol, PETSC_DECIDE, m));
 85:   PetscCall(VecSetFromOptions(init_sol));

 87:   /* create vector to hold true ts soln for comparison */
 88:   PetscCall(VecDuplicate(init_sol, &appctx.solution));

 90:   /* create LHS matrix Amat */
 91:   /*------------------------*/
 92:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, m, m, 3, NULL, &appctx.Amat));
 93:   PetscCall(MatSetFromOptions(appctx.Amat));
 94:   PetscCall(MatSetUp(appctx.Amat));
 95:   /* set space grid points - interio points only! */
 96:   PetscCall(PetscMalloc1(nz + 1, &z));
 97:   for (i = 0; i < nz; i++) z[i] = (i) * ((zFinal - zInitial) / (nz - 1));
 98:   appctx.z = z;
 99:   PetscCall(femA(&appctx, nz, z));

101:   /* create the jacobian matrix */
102:   /*----------------------------*/
103:   PetscCall(MatCreate(PETSC_COMM_WORLD, &Jmat));
104:   PetscCall(MatSetSizes(Jmat, PETSC_DECIDE, PETSC_DECIDE, m, m));
105:   PetscCall(MatSetFromOptions(Jmat));
106:   PetscCall(MatSetUp(Jmat));

108:   /* create working vectors for formulating rhs=inv(Alhs)*(Arhs*U + g) */
109:   PetscCall(VecDuplicate(init_sol, &appctx.ksp_rhs));
110:   PetscCall(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:     PetscCall(VecSetValue(init_sol, i, (PetscScalar)val, INSERT_VALUES));
117:   }
118:   PetscCall(VecAssemblyBegin(init_sol));
119:   PetscCall(VecAssemblyEnd(init_sol));

121:   /*create a time-stepping context and set the problem type */
122:   /*--------------------------------------------------------*/
123:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
124:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));

126:   /* set time-step method */
127:   PetscCall(TSSetType(ts, TSCN));

129:   /* Set optional user-defined monitoring routine */
130:   PetscCall(TSMonitorSet(ts, Monitor, &appctx, NULL));
131:   /* set the right hand side of U_t = RHSfunction(U,t) */
132:   PetscCall(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:     PetscCall(TSSetIFunction(ts, NULL, TSComputeIFunctionLinear, &appctx));
142:     PetscCall(TSSetIJacobian(ts, appctx.Amat, appctx.Amat, TSComputeIJacobianConstant, &appctx));
143:   }

145:   /* use petsc to compute the jacobian by finite differences */
146:   PetscCall(TSGetSNES(ts, &snes));
147:   PetscCall(SNESSetJacobian(snes, Jmat, Jmat, SNESComputeJacobianDefault, NULL));

149:   /* get the command line options if there are any and set them */
150:   PetscCall(TSSetFromOptions(ts));

152: #if defined(PETSC_HAVE_SUNDIALS2)
153:   {
154:     TSType    type;
155:     PetscBool sundialstype = PETSC_FALSE;
156:     PetscCall(TSGetType(ts, &type));
157:     PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &sundialstype));
158:     PetscCheck(!sundialstype || !appctx.useAlhs, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot use Alhs formulation for TSSUNDIALS type");
159:   }
160: #endif
161:   /* Sets the initial solution */
162:   PetscCall(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) 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)));
168:     PetscCall(TSSetTime(ts, ftime));
169:     PetscCall(TSSetTimeStep(ts, stepsz[k]));
170:     PetscCall(TSSetMaxSteps(ts, max_steps));
171:     PetscCall(TSSetMaxTime(ts, (k + 1) * T));
172:     PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

174:     /* loop over time steps */
175:     /*----------------------*/
176:     PetscCall(TSSolve(ts, init_sol));
177:     PetscCall(TSGetSolveTime(ts, &ftime));
178:     PetscCall(TSGetStepNumber(ts, &steps));
179:     stepsz[k + 1] = stepsz[k] * 1.5; /* change step size for the next phase */
180:   }

182:   /* free space */
183:   PetscCall(TSDestroy(&ts));
184:   PetscCall(MatDestroy(&appctx.Amat));
185:   PetscCall(MatDestroy(&Jmat));
186:   PetscCall(VecDestroy(&appctx.ksp_rhs));
187:   PetscCall(VecDestroy(&appctx.ksp_sol));
188:   PetscCall(VecDestroy(&init_sol));
189:   PetscCall(VecDestroy(&appctx.solution));
190:   PetscCall(PetscFree(z));

192:   PetscCall(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:   PetscFunctionBeginUser;
234:   /* Compute the exact solution */
235:   PetscCall(VecGetArrayWrite(appctx->solution, &u_exact));
236:   for (i = 0; i < m; i++) u_exact[i] = exact(appctx->z[i + 1], time);
237:   PetscCall(VecRestoreArrayWrite(appctx->solution, &u_exact));

239:   /* Print debugging information if desired */
240:   if (appctx->debug) {
241:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Computed solution vector at time %g\n", (double)time));
242:     PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF));
243:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Exact solution vector\n"));
244:     PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
245:   }

247:   /* Compute the 2-norm and max-norm of the error */
248:   PetscCall(VecAXPY(appctx->solution, -1.0, u));
249:   PetscCall(VecNorm(appctx->solution, NORM_2, &norm_2));

251:   norm_2 = PetscSqrtReal(h) * norm_2;
252:   PetscCall(VecNorm(appctx->solution, NORM_MAX, &norm_max));
253:   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));

255:   /*
256:      Print debugging information if desired
257:   */
258:   if (appctx->debug) {
259:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error vector\n"));
260:     PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
261:   }
262:   PetscFunctionReturn(PETSC_SUCCESS);
263: }

265: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
266:       Function to solve a linear system using KSP
267: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

269: PetscErrorCode Petsc_KSPSolve(AppCtx *obj)
270: {
271:   KSP ksp;
272:   PC  pc;

274:   PetscFunctionBeginUser;
275:   /*create the ksp context and set the operators,that is, associate the system matrix with it*/
276:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
277:   PetscCall(KSPSetOperators(ksp, obj->Amat, obj->Amat));

279:   /*get the preconditioner context, set its type and the tolerances*/
280:   PetscCall(KSPGetPC(ksp, &pc));
281:   PetscCall(PCSetType(pc, PCLU));
282:   PetscCall(KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));

284:   /*get the command line options if there are any and set them*/
285:   PetscCall(KSPSetFromOptions(ksp));

287:   /*get the linear system (ksp) solve*/
288:   PetscCall(KSPSolve(ksp, obj->ksp_rhs, obj->ksp_sol));

290:   PetscCall(KSPDestroy(&ksp));
291:   PetscFunctionReturn(PETSC_SUCCESS);
292: }

294: /***********************************************************************
295:   Function to return value of basis function or derivative of basis function.
296:  ***********************************************************************

298:         Arguments:
299:           x       = array of xpoints or nodal values
300:           xx      = point at which the basis function is to be
301:                       evaluated.
302:           il      = interval containing xx.
303:           iq      = indicates which of the two basis functions in
304:                       interval intrvl should be used
305:           nll     = array containing the endpoints of each interval.
306:           id      = If id ~= 2, the value of the basis function
307:                       is calculated; if id = 2, the value of the
308:                       derivative of the basis function is returned.
309:  ***********************************************************************/

311: PetscScalar bspl(PetscScalar *x, PetscScalar xx, PetscInt il, PetscInt iq, PetscInt nll[][2], PetscInt id)
312: {
313:   PetscScalar x1, x2, bfcn;
314:   PetscInt    i1, i2, iq1, iq2;

316:   /* Determine which basis function in interval intrvl is to be used in */
317:   iq1 = iq;
318:   if (iq1 == 0) iq2 = 1;
319:   else iq2 = 0;

321:   /*    Determine endpoint of the interval intrvl   */
322:   i1 = nll[il][iq1];
323:   i2 = nll[il][iq2];

325:   /*   Determine nodal values at the endpoints of the interval intrvl   */
326:   x1 = x[i1];
327:   x2 = x[i2];

329:   /*   Evaluate basis function   */
330:   if (id == 2) bfcn = (1.0) / (x1 - x2);
331:   else bfcn = (xx - x2) / (x1 - x2);
332:   return bfcn;
333: }

335: /*---------------------------------------------------------
336:   Function called by rhs function to get B and g
337: ---------------------------------------------------------*/
338: PetscErrorCode femBg(PetscScalar btri[][3], PetscScalar *f, PetscInt nz, PetscScalar *z, PetscReal t)
339: {
340:   PetscInt    i, j, jj, il, ip, ipp, ipq, iq, iquad, iqq;
341:   PetscInt    nli[num_z][2], indx[num_z];
342:   PetscScalar dd, dl, zip, zipq, zz, b_z, bb_z, bij;
343:   PetscScalar zquad[num_z][3], dlen[num_z], qdwt[3];

345:   PetscFunctionBeginUser;
346:   /*  initializing everything - btri and f are initialized in rhs.c  */
347:   for (i = 0; i < nz; i++) {
348:     nli[i][0]   = 0;
349:     nli[i][1]   = 0;
350:     indx[i]     = 0;
351:     zquad[i][0] = 0.0;
352:     zquad[i][1] = 0.0;
353:     zquad[i][2] = 0.0;
354:     dlen[i]     = 0.0;
355:   } /*end for (i)*/

357:   /*  quadrature weights  */
358:   qdwt[0] = 1.0 / 6.0;
359:   qdwt[1] = 4.0 / 6.0;
360:   qdwt[2] = 1.0 / 6.0;

362:   /* 1st and last nodes have Dirichlet boundary condition -
363:      set indices there to -1 */

365:   for (i = 0; i < nz - 1; i++) indx[i] = i - 1;
366:   indx[nz - 1] = -1;

368:   ipq = 0;
369:   for (il = 0; il < nz - 1; il++) {
370:     ip           = ipq;
371:     ipq          = ip + 1;
372:     zip          = z[ip];
373:     zipq         = z[ipq];
374:     dl           = zipq - zip;
375:     zquad[il][0] = zip;
376:     zquad[il][1] = (0.5) * (zip + zipq);
377:     zquad[il][2] = zipq;
378:     dlen[il]     = PetscAbsScalar(dl);
379:     nli[il][0]   = ip;
380:     nli[il][1]   = ipq;
381:   }

383:   for (il = 0; il < nz - 1; il++) {
384:     for (iquad = 0; iquad < 3; iquad++) {
385:       dd = (dlen[il]) * (qdwt[iquad]);
386:       zz = zquad[il][iquad];

388:       for (iq = 0; iq < 2; iq++) {
389:         ip  = nli[il][iq];
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:             bb_z = bspl(z, zz, il, iqq, nli, 2);
397:             j    = indx[ipp];
398:             bij  = -b_z * bb_z;

400:             if (j > -1) {
401:               jj = 1 + j - i;
402:               btri[i][jj] += bij * dd;
403:             } else {
404:               f[i] += bij * dd * exact(z[ipp], t);
405:               /* f[i] += 0.0; */
406:               /* if (il==0 && j==-1) { */
407:               /* f[i] += bij*dd*exact(zz,t); */
408:               /* }*/ /*end if*/
409:             }        /*end else*/
410:           }          /*end for (iqq)*/
411:         }            /*end if (i>0)*/
412:       }              /*end for (iq)*/
413:     }                /*end for (iquad)*/
414:   }                  /*end for (il)*/
415:   PetscFunctionReturn(PETSC_SUCCESS);
416: }

418: PetscErrorCode femA(AppCtx *obj, PetscInt nz, PetscScalar *z)
419: {
420:   PetscInt    i, j, il, ip, ipp, ipq, iq, iquad, iqq;
421:   PetscInt    nli[num_z][2], indx[num_z];
422:   PetscScalar dd, dl, zip, zipq, zz, bb, bbb, aij;
423:   PetscScalar rquad[num_z][3], dlen[num_z], qdwt[3], add_term;

425:   PetscFunctionBeginUser;
426:   /*  initializing everything  */
427:   for (i = 0; i < nz; i++) {
428:     nli[i][0]   = 0;
429:     nli[i][1]   = 0;
430:     indx[i]     = 0;
431:     rquad[i][0] = 0.0;
432:     rquad[i][1] = 0.0;
433:     rquad[i][2] = 0.0;
434:     dlen[i]     = 0.0;
435:   } /*end for (i)*/

437:   /*  quadrature weights  */
438:   qdwt[0] = 1.0 / 6.0;
439:   qdwt[1] = 4.0 / 6.0;
440:   qdwt[2] = 1.0 / 6.0;

442:   /* 1st and last nodes have Dirichlet boundary condition -
443:      set indices there to -1 */

445:   for (i = 0; i < nz - 1; i++) indx[i] = i - 1;
446:   indx[nz - 1] = -1;

448:   ipq = 0;

450:   for (il = 0; il < nz - 1; il++) {
451:     ip           = ipq;
452:     ipq          = ip + 1;
453:     zip          = z[ip];
454:     zipq         = z[ipq];
455:     dl           = zipq - zip;
456:     rquad[il][0] = zip;
457:     rquad[il][1] = (0.5) * (zip + zipq);
458:     rquad[il][2] = zipq;
459:     dlen[il]     = PetscAbsScalar(dl);
460:     nli[il][0]   = ip;
461:     nli[il][1]   = ipq;
462:   } /*end for (il)*/

464:   for (il = 0; il < nz - 1; il++) {
465:     for (iquad = 0; iquad < 3; iquad++) {
466:       dd = (dlen[il]) * (qdwt[iquad]);
467:       zz = rquad[il][iquad];

469:       for (iq = 0; iq < 2; iq++) {
470:         ip = nli[il][iq];
471:         bb = bspl(z, zz, il, iq, nli, 1);
472:         i  = indx[ip];
473:         if (i > -1) {
474:           for (iqq = 0; iqq < 2; iqq++) {
475:             ipp = nli[il][iqq];
476:             bbb = bspl(z, zz, il, iqq, nli, 1);
477:             j   = indx[ipp];
478:             aij = bb * bbb;
479:             if (j > -1) {
480:               add_term = aij * dd;
481:               PetscCall(MatSetValue(obj->Amat, i, j, add_term, ADD_VALUES));
482:             } /*endif*/
483:           }   /*end for (iqq)*/
484:         }     /*end if (i>0)*/
485:       }       /*end for (iq)*/
486:     }         /*end for (iquad)*/
487:   }           /*end for (il)*/
488:   PetscCall(MatAssemblyBegin(obj->Amat, MAT_FINAL_ASSEMBLY));
489:   PetscCall(MatAssemblyEnd(obj->Amat, MAT_FINAL_ASSEMBLY));
490:   PetscFunctionReturn(PETSC_SUCCESS);
491: }

493: /*---------------------------------------------------------
494:         Function to fill the rhs vector with
495:         By + g values ****
496: ---------------------------------------------------------*/
497: PetscErrorCode rhs(AppCtx *obj, PetscScalar *y, PetscInt nz, PetscScalar *z, PetscReal t)
498: {
499:   PetscInt    i, j, js, je, jj;
500:   PetscScalar val, g[num_z], btri[num_z][3], add_term;

502:   PetscFunctionBeginUser;
503:   for (i = 0; i < nz - 2; i++) {
504:     for (j = 0; j <= 2; j++) btri[i][j] = 0.0;
505:     g[i] = 0.0;
506:   }

508:   /*  call femBg to set the tri-diagonal b matrix and vector g  */
509:   PetscCall(femBg(btri, g, nz, z, t));

511:   /*  setting the entries of the right hand side vector  */
512:   for (i = 0; i < nz - 2; i++) {
513:     val = 0.0;
514:     js  = 0;
515:     if (i == 0) js = 1;
516:     je = 2;
517:     if (i == nz - 2) je = 1;

519:     for (jj = js; jj <= je; jj++) {
520:       j = i + jj - 1;
521:       val += (btri[i][jj]) * (y[j]);
522:     }
523:     add_term = val + g[i];
524:     PetscCall(VecSetValue(obj->ksp_rhs, (PetscInt)i, (PetscScalar)add_term, INSERT_VALUES));
525:   }
526:   PetscCall(VecAssemblyBegin(obj->ksp_rhs));
527:   PetscCall(VecAssemblyEnd(obj->ksp_rhs));
528:   PetscFunctionReturn(PETSC_SUCCESS);
529: }

531: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
532: %%   Function to form the right hand side of the time-stepping problem.                       %%
533: %% -------------------------------------------------------------------------------------------%%
534:   if (useAlhs):
535:     globalout = By+g
536:   else if (!useAlhs):
537:     globalout = f(y,t)=Ainv(By+g),
538:       in which the ksp solver to transform the problem A*ydot=By+g
539:       to the problem ydot=f(y,t)=inv(A)*(By+g)
540: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

542: PetscErrorCode RHSfunction(TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx)
543: {
544:   AppCtx            *obj = (AppCtx *)ctx;
545:   PetscScalar        soln[num_z];
546:   const PetscScalar *soln_ptr;
547:   PetscInt           i, nz = obj->nz;
548:   PetscReal          time;

550:   PetscFunctionBeginUser;
551:   /* get the previous solution to compute updated system */
552:   PetscCall(VecGetArrayRead(globalin, &soln_ptr));
553:   for (i = 0; i < num_z - 2; i++) soln[i] = soln_ptr[i];
554:   PetscCall(VecRestoreArrayRead(globalin, &soln_ptr));
555:   soln[num_z - 1] = 0.0;
556:   soln[num_z - 2] = 0.0;

558:   /* clear out the matrix and rhs for ksp to keep things straight */
559:   PetscCall(VecSet(obj->ksp_rhs, (PetscScalar)0.0));

561:   time = t;
562:   /* get the updated system */
563:   PetscCall(rhs(obj, soln, nz, obj->z, time)); /* setup of the By+g rhs */

565:   /* do a ksp solve to get the rhs for the ts problem */
566:   if (obj->useAlhs) {
567:     /* ksp_sol = ksp_rhs */
568:     PetscCall(VecCopy(obj->ksp_rhs, globalout));
569:   } else {
570:     /* ksp_sol = inv(Amat)*ksp_rhs */
571:     PetscCall(Petsc_KSPSolve(obj));
572:     PetscCall(VecCopy(obj->ksp_sol, globalout));
573:   }
574:   PetscFunctionReturn(PETSC_SUCCESS);
575: }

577: /*TEST

579:     build:
580:       requires: !complex

582:     test:
583:       suffix: euler
584:       output_file: output/ex3.out

586:     test:
587:       suffix: 2
588:       args:   -useAlhs
589:       output_file: output/ex3.out
590:       TODO: Broken because SNESComputeJacobianDefault is incompatible with TSComputeIJacobianConstant

592: TEST*/