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";
5
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 ----------------------------------------------------------------------------*/
11
12 #include <petscksp.h>
13 #include <petscts.h>
14
15 /* special variable - max size of all arrays */
16 #define num_z 10
17
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;
33
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 *);
42
main(int argc,char ** argv)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;
54
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");
59
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 = 10000;
67
68 appctx.m = m;
69 appctx.max_probsz = nz;
70 appctx.debug = PETSC_FALSE;
71 appctx.useAlhs = PETSC_FALSE;
72
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;
79
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));
85
86 /* create vector to hold true ts soln for comparison */
87 PetscCall(VecDuplicate(init_sol, &appctx.solution));
88
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));
99
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));
106
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));
110
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, val, INSERT_VALUES));
116 }
117 PetscCall(VecAssemblyBegin(init_sol));
118 PetscCall(VecAssemblyEnd(init_sol));
119
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));
124
125 /* set time-step method */
126 PetscCall(TSSetType(ts, TSCN));
127
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));
132
133 if (appctx.useAlhs) {
134 /* set the left hand side matrix of Amat*U_t = rhs(U,t) */
135
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 }
143
144 /* use PETSc to compute the jacobian by finite differences */
145 PetscCall(TSGetSNES(ts, &snes));
146 PetscCall(SNESSetJacobian(snes, Jmat, Jmat, SNESComputeJacobianDefault, NULL));
147
148 /* get the command line options if there are any and set them */
149 PetscCall(TSSetFromOptions(ts));
150
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));
162
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));
172
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 }
180
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));
190
191 PetscCall(PetscFinalize());
192 return 0;
193 }
194
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 --------------------------------------------------------------------------*/
exact(PetscScalar z,PetscReal t)199 PetscScalar exact(PetscScalar z, PetscReal t)
200 {
201 PetscScalar val, ex1, ex2;
202
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 }
208
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.
213
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 */
Monitor(TS ts,PetscInt step,PetscReal time,Vec u,PetscCtx ctx)225 PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec u, PetscCtx 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;
231
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));
237
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 }
245
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));
249
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));
253
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 }
263
264 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
265 Function to solve a linear system using KSP
266 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
267
Petsc_KSPSolve(AppCtx * obj)268 PetscErrorCode Petsc_KSPSolve(AppCtx *obj)
269 {
270 KSP ksp;
271 PC pc;
272
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));
277
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));
282
283 /*get the command line options if there are any and set them*/
284 PetscCall(KSPSetFromOptions(ksp));
285
286 /*get the linear system (ksp) solve*/
287 PetscCall(KSPSolve(ksp, obj->ksp_rhs, obj->ksp_sol));
288
289 PetscCall(KSPDestroy(&ksp));
290 PetscFunctionReturn(PETSC_SUCCESS);
291 }
292
293 /***********************************************************************
294 Function to return value of basis function or derivative of basis function.
295 ***********************************************************************
296
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 ***********************************************************************/
309
bspl(PetscScalar * x,PetscScalar xx,PetscInt il,PetscInt iq,PetscInt nll[][2],PetscInt id)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;
314
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;
319
320 /* Determine endpoint of the interval intrvl */
321 i1 = nll[il][iq1];
322 i2 = nll[il][iq2];
323
324 /* Determine nodal values at the endpoints of the interval intrvl */
325 x1 = x[i1];
326 x2 = x[i2];
327
328 /* Evaluate basis function */
329 if (id == 2) bfcn = (1.0) / (x1 - x2);
330 else bfcn = (xx - x2) / (x1 - x2);
331 return bfcn;
332 }
333
334 /*---------------------------------------------------------
335 Function called by rhs function to get B and g
336 ---------------------------------------------------------*/
femBg(PetscScalar btri[][3],PetscScalar * f,PetscInt nz,PetscScalar * z,PetscReal t)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];
343
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)*/
355
356 /* quadrature weights */
357 qdwt[0] = 1.0 / 6.0;
358 qdwt[1] = 4.0 / 6.0;
359 qdwt[2] = 1.0 / 6.0;
360
361 /* 1st and last nodes have Dirichlet boundary condition -
362 set indices there to -1 */
363
364 for (i = 0; i < nz - 1; i++) indx[i] = i - 1;
365 indx[nz - 1] = -1;
366
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 }
381
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];
386
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];
391
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;
398
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 }
416
femA(AppCtx * obj,PetscInt nz,PetscScalar * z)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;
423
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)*/
435
436 /* quadrature weights */
437 qdwt[0] = 1.0 / 6.0;
438 qdwt[1] = 4.0 / 6.0;
439 qdwt[2] = 1.0 / 6.0;
440
441 /* 1st and last nodes have Dirichlet boundary condition -
442 set indices there to -1 */
443
444 for (i = 0; i < nz - 1; i++) indx[i] = i - 1;
445 indx[nz - 1] = -1;
446
447 ipq = 0;
448
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)*/
462
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];
467
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 }
491
492 /*---------------------------------------------------------
493 Function to fill the rhs vector with
494 By + g values ****
495 ---------------------------------------------------------*/
rhs(AppCtx * obj,PetscScalar * y,PetscInt nz,PetscScalar * z,PetscReal t)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;
500
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 }
506
507 /* call femBg to set the tri-diagonal b matrix and vector g */
508 PetscCall(femBg(btri, g, nz, z, t));
509
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;
517
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, i, add_term, INSERT_VALUES));
524 }
525 PetscCall(VecAssemblyBegin(obj->ksp_rhs));
526 PetscCall(VecAssemblyEnd(obj->ksp_rhs));
527 PetscFunctionReturn(PETSC_SUCCESS);
528 }
529
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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
540
RHSfunction(TS ts,PetscReal t,Vec globalin,Vec globalout,PetscCtx ctx)541 PetscErrorCode RHSfunction(TS ts, PetscReal t, Vec globalin, Vec globalout, PetscCtx 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;
548
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;
556
557 /* clear out the matrix and rhs for ksp to keep things straight */
558 PetscCall(VecSet(obj->ksp_rhs, (PetscScalar)0.0));
559
560 time = t;
561 /* get the updated system */
562 PetscCall(rhs(obj, soln, nz, obj->z, time)); /* setup of the By+g rhs */
563
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 }
575
576 /*TEST
577
578 build:
579 requires: !complex
580
581 test:
582 suffix: euler
583 output_file: output/ex3.out
584
585 test:
586 suffix: 2
587 args: -useAlhs
588 output_file: output/ex3.out
589 TODO: Broken because SNESComputeJacobianDefault is incompatible with TSComputeIJacobianConstant
590
591 TEST*/
592