1
2 static char help[] ="Solves a simple time-dependent linear PDE (the heat equation).\n\
3 Input parameters include:\n\
4 -m <points>, where <points> = number of grid points\n\
5 -time_dependent_rhs : Treat the problem as having a time-dependent right-hand side\n\
6 -use_ifunc : Use IFunction/IJacobian interface\n\
7 -debug : Activate debugging printouts\n\
8 -nox : Deactivate x-window graphics\n\n";
9
10 /*
11 Concepts: TS^time-dependent linear problems
12 Concepts: TS^heat equation
13 Concepts: TS^diffusion equation
14 Processors: 1
15 */
16
17 /* ------------------------------------------------------------------------
18
19 This program solves the one-dimensional heat equation (also called the
20 diffusion equation),
21 u_t = u_xx,
22 on the domain 0 <= x <= 1, with the boundary conditions
23 u(t,0) = 0, u(t,1) = 0,
24 and the initial condition
25 u(0,x) = sin(6*pi*x) + 3*sin(2*pi*x).
26 This is a linear, second-order, parabolic equation.
27
28 We discretize the right-hand side using finite differences with
29 uniform grid spacing h:
30 u_xx = (u_{i+1} - 2u_{i} + u_{i-1})/(h^2)
31 We then demonstrate time evolution using the various TS methods by
32 running the program via
33 ex3 -ts_type <timestepping solver>
34
35 We compare the approximate solution with the exact solution, given by
36 u_exact(x,t) = exp(-36*pi*pi*t) * sin(6*pi*x) +
37 3*exp(-4*pi*pi*t) * sin(2*pi*x)
38
39 Notes:
40 This code demonstrates the TS solver interface to two variants of
41 linear problems, u_t = f(u,t), namely
42 - time-dependent f: f(u,t) is a function of t
43 - time-independent f: f(u,t) is simply f(u)
44
45 The parallel version of this code is ts/tutorials/ex4.c
46
47 ------------------------------------------------------------------------- */
48
49 /*
50 Include "petscts.h" so that we can use TS solvers. Note that this file
51 automatically includes:
52 petscsys.h - base PETSc routines petscvec.h - vectors
53 petscmat.h - matrices
54 petscis.h - index sets petscksp.h - Krylov subspace methods
55 petscviewer.h - viewers petscpc.h - preconditioners
56 petscksp.h - linear solvers petscsnes.h - nonlinear solvers
57 */
58
59 #include <petscts.h>
60 #include <petscdraw.h>
61
62 /*
63 User-defined application context - contains data needed by the
64 application-provided call-back routines.
65 */
66 typedef struct {
67 Vec solution; /* global exact solution vector */
68 PetscInt m; /* total number of grid points */
69 PetscReal h; /* mesh width h = 1/(m-1) */
70 PetscBool debug; /* flag (1 indicates activation of debugging printouts) */
71 PetscViewer viewer1,viewer2; /* viewers for the solution and error */
72 PetscReal norm_2,norm_max; /* error norms */
73 Mat A; /* RHS mat, used with IFunction interface */
74 PetscReal oshift; /* old shift applied, prevent to recompute the IJacobian */
75 } AppCtx;
76
77 /*
78 User-defined routines
79 */
80 extern PetscErrorCode InitialConditions(Vec,AppCtx*);
81 extern PetscErrorCode RHSMatrixHeat(TS,PetscReal,Vec,Mat,Mat,void*);
82 extern PetscErrorCode IFunctionHeat(TS,PetscReal,Vec,Vec,Vec,void*);
83 extern PetscErrorCode IJacobianHeat(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
84 extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);
85 extern PetscErrorCode ExactSolution(PetscReal,Vec,AppCtx*);
86
main(int argc,char ** argv)87 int main(int argc,char **argv)
88 {
89 AppCtx appctx; /* user-defined application context */
90 TS ts; /* timestepping context */
91 Mat A; /* matrix data structure */
92 Vec u; /* approximate solution vector */
93 PetscReal time_total_max = 100.0; /* default max total time */
94 PetscInt time_steps_max = 100; /* default max timesteps */
95 PetscDraw draw; /* drawing context */
96 PetscErrorCode ierr;
97 PetscInt steps,m;
98 PetscMPIInt size;
99 PetscReal dt;
100 PetscBool flg,flg_string;
101
102 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103 Initialize program and set problem parameters
104 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105
106 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
107 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
108 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");
109
110 m = 60;
111 ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr);
112 ierr = PetscOptionsHasName(NULL,NULL,"-debug",&appctx.debug);CHKERRQ(ierr);
113 flg_string = PETSC_FALSE;
114 ierr = PetscOptionsGetBool(NULL,NULL,"-test_string_viewer",&flg_string,NULL);CHKERRQ(ierr);
115
116 appctx.m = m;
117 appctx.h = 1.0/(m-1.0);
118 appctx.norm_2 = 0.0;
119 appctx.norm_max = 0.0;
120
121 ierr = PetscPrintf(PETSC_COMM_SELF,"Solving a linear TS problem on 1 processor\n");CHKERRQ(ierr);
122
123 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124 Create vector data structures
125 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126
127 /*
128 Create vector data structures for approximate and exact solutions
129 */
130 ierr = VecCreateSeq(PETSC_COMM_SELF,m,&u);CHKERRQ(ierr);
131 ierr = VecDuplicate(u,&appctx.solution);CHKERRQ(ierr);
132
133 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134 Set up displays to show graphs of the solution and error
135 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136
137 ierr = PetscViewerDrawOpen(PETSC_COMM_SELF,0,"",80,380,400,160,&appctx.viewer1);CHKERRQ(ierr);
138 ierr = PetscViewerDrawGetDraw(appctx.viewer1,0,&draw);CHKERRQ(ierr);
139 ierr = PetscDrawSetDoubleBuffer(draw);CHKERRQ(ierr);
140 ierr = PetscViewerDrawOpen(PETSC_COMM_SELF,0,"",80,0,400,160,&appctx.viewer2);CHKERRQ(ierr);
141 ierr = PetscViewerDrawGetDraw(appctx.viewer2,0,&draw);CHKERRQ(ierr);
142 ierr = PetscDrawSetDoubleBuffer(draw);CHKERRQ(ierr);
143
144 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145 Create timestepping solver context
146 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147
148 ierr = TSCreate(PETSC_COMM_SELF,&ts);CHKERRQ(ierr);
149 ierr = TSSetProblemType(ts,TS_LINEAR);CHKERRQ(ierr);
150
151 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152 Set optional user-defined monitoring routine
153 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154
155 if (!flg_string) {
156 ierr = TSMonitorSet(ts,Monitor,&appctx,NULL);CHKERRQ(ierr);
157 }
158
159 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160
161 Create matrix data structure; set matrix evaluation routine.
162 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
163
164 ierr = MatCreate(PETSC_COMM_SELF,&A);CHKERRQ(ierr);
165 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,m);CHKERRQ(ierr);
166 ierr = MatSetFromOptions(A);CHKERRQ(ierr);
167 ierr = MatSetUp(A);CHKERRQ(ierr);
168
169 flg = PETSC_FALSE;
170 ierr = PetscOptionsGetBool(NULL,NULL,"-use_ifunc",&flg,NULL);CHKERRQ(ierr);
171 if (!flg) {
172 appctx.A = NULL;
173 ierr = PetscOptionsGetBool(NULL,NULL,"-time_dependent_rhs",&flg,NULL);CHKERRQ(ierr);
174 if (flg) {
175 /*
176 For linear problems with a time-dependent f(u,t) in the equation
177 u_t = f(u,t), the user provides the discretized right-hand-side
178 as a time-dependent matrix.
179 */
180 ierr = TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx);CHKERRQ(ierr);
181 ierr = TSSetRHSJacobian(ts,A,A,RHSMatrixHeat,&appctx);CHKERRQ(ierr);
182 } else {
183 /*
184 For linear problems with a time-independent f(u) in the equation
185 u_t = f(u), the user provides the discretized right-hand-side
186 as a matrix only once, and then sets the special Jacobian evaluation
187 routine TSComputeRHSJacobianConstant() which will NOT recompute the Jacobian.
188 */
189 ierr = RHSMatrixHeat(ts,0.0,u,A,A,&appctx);CHKERRQ(ierr);
190 ierr = TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx);CHKERRQ(ierr);
191 ierr = TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,&appctx);CHKERRQ(ierr);
192 }
193 } else {
194 Mat J;
195
196 ierr = RHSMatrixHeat(ts,0.0,u,A,A,&appctx);CHKERRQ(ierr);
197 ierr = MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&J);CHKERRQ(ierr);
198 ierr = TSSetIFunction(ts,NULL,IFunctionHeat,&appctx);CHKERRQ(ierr);
199 ierr = TSSetIJacobian(ts,J,J,IJacobianHeat,&appctx);CHKERRQ(ierr);
200 ierr = MatDestroy(&J);CHKERRQ(ierr);
201
202 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
203 appctx.A = A;
204 appctx.oshift = PETSC_MIN_REAL;
205 }
206 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
207 Set solution vector and initial timestep
208 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
209
210 dt = appctx.h*appctx.h/2.0;
211 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
212
213 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
214 Customize timestepping solver:
215 - Set the solution method to be the Backward Euler method.
216 - Set timestepping duration info
217 Then set runtime options, which can override these defaults.
218 For example,
219 -ts_max_steps <maxsteps> -ts_max_time <maxtime>
220 to override the defaults set by TSSetMaxSteps()/TSSetMaxTime().
221 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
222
223 ierr = TSSetMaxSteps(ts,time_steps_max);CHKERRQ(ierr);
224 ierr = TSSetMaxTime(ts,time_total_max);CHKERRQ(ierr);
225 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
226 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
227
228 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
229 Solve the problem
230 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
231
232 /*
233 Evaluate initial conditions
234 */
235 ierr = InitialConditions(u,&appctx);CHKERRQ(ierr);
236
237 /*
238 Run the timestepping solver
239 */
240 ierr = TSSolve(ts,u);CHKERRQ(ierr);
241 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr);
242
243 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
244 View timestepping solver info
245 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
246
247 ierr = PetscPrintf(PETSC_COMM_SELF,"avg. error (2 norm) = %g, avg. error (max norm) = %g\n",(double)(appctx.norm_2/steps),(double)(appctx.norm_max/steps));CHKERRQ(ierr);
248 if (!flg_string) {
249 ierr = TSView(ts,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
250 } else {
251 PetscViewer stringviewer;
252 char string[512];
253 const char *outstring;
254
255 ierr = PetscViewerStringOpen(PETSC_COMM_WORLD,string,sizeof(string),&stringviewer);CHKERRQ(ierr);
256 ierr = TSView(ts,stringviewer);CHKERRQ(ierr);
257 ierr = PetscViewerStringGetStringRead(stringviewer,&outstring,NULL);CHKERRQ(ierr);
258 if ((char*)outstring != (char*)string) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"String returned from viewer does not equal original string");
259 ierr = PetscPrintf(PETSC_COMM_WORLD,"Output from string viewer:%s\n",outstring);CHKERRQ(ierr);
260 ierr = PetscViewerDestroy(&stringviewer);CHKERRQ(ierr);
261 }
262
263 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
264 Free work space. All PETSc objects should be destroyed when they
265 are no longer needed.
266 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
267
268 ierr = TSDestroy(&ts);CHKERRQ(ierr);
269 ierr = MatDestroy(&A);CHKERRQ(ierr);
270 ierr = VecDestroy(&u);CHKERRQ(ierr);
271 ierr = PetscViewerDestroy(&appctx.viewer1);CHKERRQ(ierr);
272 ierr = PetscViewerDestroy(&appctx.viewer2);CHKERRQ(ierr);
273 ierr = VecDestroy(&appctx.solution);CHKERRQ(ierr);
274 ierr = MatDestroy(&appctx.A);CHKERRQ(ierr);
275
276 /*
277 Always call PetscFinalize() before exiting a program. This routine
278 - finalizes the PETSc libraries as well as MPI
279 - provides summary and diagnostic information if certain runtime
280 options are chosen (e.g., -log_view).
281 */
282 ierr = PetscFinalize();
283 return ierr;
284 }
285 /* --------------------------------------------------------------------- */
286 /*
287 InitialConditions - Computes the solution at the initial time.
288
289 Input Parameter:
290 u - uninitialized solution vector (global)
291 appctx - user-defined application context
292
293 Output Parameter:
294 u - vector with solution at initial time (global)
295 */
InitialConditions(Vec u,AppCtx * appctx)296 PetscErrorCode InitialConditions(Vec u,AppCtx *appctx)
297 {
298 PetscScalar *u_localptr,h = appctx->h;
299 PetscErrorCode ierr;
300 PetscInt i;
301
302 /*
303 Get a pointer to vector data.
304 - For default PETSc vectors, VecGetArray() returns a pointer to
305 the data array. Otherwise, the routine is implementation dependent.
306 - You MUST call VecRestoreArray() when you no longer need access to
307 the array.
308 - Note that the Fortran interface to VecGetArray() differs from the
309 C version. See the users manual for details.
310 */
311 ierr = VecGetArrayWrite(u,&u_localptr);CHKERRQ(ierr);
312
313 /*
314 We initialize the solution array by simply writing the solution
315 directly into the array locations. Alternatively, we could use
316 VecSetValues() or VecSetValuesLocal().
317 */
318 for (i=0; i<appctx->m; i++) u_localptr[i] = PetscSinScalar(PETSC_PI*i*6.*h) + 3.*PetscSinScalar(PETSC_PI*i*2.*h);
319
320 /*
321 Restore vector
322 */
323 ierr = VecRestoreArrayWrite(u,&u_localptr);CHKERRQ(ierr);
324
325 /*
326 Print debugging information if desired
327 */
328 if (appctx->debug) {
329 ierr = PetscPrintf(PETSC_COMM_WORLD,"Initial guess vector\n");CHKERRQ(ierr);
330 ierr = VecView(u,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
331 }
332
333 return 0;
334 }
335 /* --------------------------------------------------------------------- */
336 /*
337 ExactSolution - Computes the exact solution at a given time.
338
339 Input Parameters:
340 t - current time
341 solution - vector in which exact solution will be computed
342 appctx - user-defined application context
343
344 Output Parameter:
345 solution - vector with the newly computed exact solution
346 */
ExactSolution(PetscReal t,Vec solution,AppCtx * appctx)347 PetscErrorCode ExactSolution(PetscReal t,Vec solution,AppCtx *appctx)
348 {
349 PetscScalar *s_localptr,h = appctx->h,ex1,ex2,sc1,sc2,tc = t;
350 PetscErrorCode ierr;
351 PetscInt i;
352
353 /*
354 Get a pointer to vector data.
355 */
356 ierr = VecGetArrayWrite(solution,&s_localptr);CHKERRQ(ierr);
357
358 /*
359 Simply write the solution directly into the array locations.
360 Alternatively, we culd use VecSetValues() or VecSetValuesLocal().
361 */
362 ex1 = PetscExpScalar(-36.*PETSC_PI*PETSC_PI*tc);
363 ex2 = PetscExpScalar(-4.*PETSC_PI*PETSC_PI*tc);
364 sc1 = PETSC_PI*6.*h; sc2 = PETSC_PI*2.*h;
365 for (i=0; i<appctx->m; i++) s_localptr[i] = PetscSinScalar(sc1*(PetscReal)i)*ex1 + 3.*PetscSinScalar(sc2*(PetscReal)i)*ex2;
366
367 /*
368 Restore vector
369 */
370 ierr = VecRestoreArrayWrite(solution,&s_localptr);CHKERRQ(ierr);
371 return 0;
372 }
373 /* --------------------------------------------------------------------- */
374 /*
375 Monitor - User-provided routine to monitor the solution computed at
376 each timestep. This example plots the solution and computes the
377 error in two different norms.
378
379 This example also demonstrates changing the timestep via TSSetTimeStep().
380
381 Input Parameters:
382 ts - the timestep context
383 step - the count of the current step (with 0 meaning the
384 initial condition)
385 time - the current time
386 u - the solution at this timestep
387 ctx - the user-provided context for this monitoring routine.
388 In this case we use the application context which contains
389 information about the problem size, workspace and the exact
390 solution.
391 */
Monitor(TS ts,PetscInt step,PetscReal time,Vec u,void * ctx)392 PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec u,void *ctx)
393 {
394 AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */
395 PetscErrorCode ierr;
396 PetscReal norm_2,norm_max,dt,dttol;
397
398 /*
399 View a graph of the current iterate
400 */
401 ierr = VecView(u,appctx->viewer2);CHKERRQ(ierr);
402
403 /*
404 Compute the exact solution
405 */
406 ierr = ExactSolution(time,appctx->solution,appctx);CHKERRQ(ierr);
407
408 /*
409 Print debugging information if desired
410 */
411 if (appctx->debug) {
412 ierr = PetscPrintf(PETSC_COMM_SELF,"Computed solution vector\n");CHKERRQ(ierr);
413 ierr = VecView(u,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
414 ierr = PetscPrintf(PETSC_COMM_SELF,"Exact solution vector\n");CHKERRQ(ierr);
415 ierr = VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
416 }
417
418 /*
419 Compute the 2-norm and max-norm of the error
420 */
421 ierr = VecAXPY(appctx->solution,-1.0,u);CHKERRQ(ierr);
422 ierr = VecNorm(appctx->solution,NORM_2,&norm_2);CHKERRQ(ierr);
423 norm_2 = PetscSqrtReal(appctx->h)*norm_2;
424 ierr = VecNorm(appctx->solution,NORM_MAX,&norm_max);CHKERRQ(ierr);
425
426 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
427 ierr = PetscPrintf(PETSC_COMM_WORLD,"Timestep %3D: step size = %g, time = %g, 2-norm error = %g, max norm error = %g\n",step,(double)dt,(double)time,(double)norm_2,(double)norm_max);CHKERRQ(ierr);
428
429 appctx->norm_2 += norm_2;
430 appctx->norm_max += norm_max;
431
432 dttol = .0001;
433 ierr = PetscOptionsGetReal(NULL,NULL,"-dttol",&dttol,NULL);CHKERRQ(ierr);
434 if (dt < dttol) {
435 dt *= .999;
436 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
437 }
438
439 /*
440 View a graph of the error
441 */
442 ierr = VecView(appctx->solution,appctx->viewer1);CHKERRQ(ierr);
443
444 /*
445 Print debugging information if desired
446 */
447 if (appctx->debug) {
448 ierr = PetscPrintf(PETSC_COMM_SELF,"Error vector\n");CHKERRQ(ierr);
449 ierr = VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
450 }
451
452 return 0;
453 }
454 /* --------------------------------------------------------------------- */
455 /*
456 RHSMatrixHeat - User-provided routine to compute the right-hand-side
457 matrix for the heat equation.
458
459 Input Parameters:
460 ts - the TS context
461 t - current time
462 global_in - global input vector
463 dummy - optional user-defined context, as set by TSetRHSJacobian()
464
465 Output Parameters:
466 AA - Jacobian matrix
467 BB - optionally different preconditioning matrix
468 str - flag indicating matrix structure
469
470 Notes:
471 Recall that MatSetValues() uses 0-based row and column numbers
472 in Fortran as well as in C.
473 */
RHSMatrixHeat(TS ts,PetscReal t,Vec X,Mat AA,Mat BB,void * ctx)474 PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Vec X,Mat AA,Mat BB,void *ctx)
475 {
476 Mat A = AA; /* Jacobian matrix */
477 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
478 PetscInt mstart = 0;
479 PetscInt mend = appctx->m;
480 PetscErrorCode ierr;
481 PetscInt i,idx[3];
482 PetscScalar v[3],stwo = -2./(appctx->h*appctx->h),sone = -.5*stwo;
483
484 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
485 Compute entries for the locally owned part of the matrix
486 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
487 /*
488 Set matrix rows corresponding to boundary data
489 */
490
491 mstart = 0;
492 v[0] = 1.0;
493 ierr = MatSetValues(A,1,&mstart,1,&mstart,v,INSERT_VALUES);CHKERRQ(ierr);
494 mstart++;
495
496 mend--;
497 v[0] = 1.0;
498 ierr = MatSetValues(A,1,&mend,1,&mend,v,INSERT_VALUES);CHKERRQ(ierr);
499
500 /*
501 Set matrix rows corresponding to interior data. We construct the
502 matrix one row at a time.
503 */
504 v[0] = sone; v[1] = stwo; v[2] = sone;
505 for (i=mstart; i<mend; i++) {
506 idx[0] = i-1; idx[1] = i; idx[2] = i+1;
507 ierr = MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);CHKERRQ(ierr);
508 }
509
510 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
511 Complete the matrix assembly process and set some options
512 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
513 /*
514 Assemble matrix, using the 2-step process:
515 MatAssemblyBegin(), MatAssemblyEnd()
516 Computations can be done while messages are in transition
517 by placing code between these two statements.
518 */
519 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
520 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
521
522 /*
523 Set and option to indicate that we will never add a new nonzero location
524 to the matrix. If we do, it will generate an error.
525 */
526 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
527
528 return 0;
529 }
530
IFunctionHeat(TS ts,PetscReal t,Vec X,Vec Xdot,Vec r,void * ctx)531 PetscErrorCode IFunctionHeat(TS ts,PetscReal t,Vec X,Vec Xdot,Vec r,void *ctx)
532 {
533 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
534 PetscErrorCode ierr;
535
536 ierr = MatMult(appctx->A,X,r);CHKERRQ(ierr);
537 ierr = VecAYPX(r,-1.0,Xdot);CHKERRQ(ierr);
538 return 0;
539 }
540
IJacobianHeat(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal s,Mat A,Mat B,void * ctx)541 PetscErrorCode IJacobianHeat(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal s,Mat A,Mat B,void *ctx)
542 {
543 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
544 PetscErrorCode ierr;
545
546 if (appctx->oshift == s) return 0;
547 ierr = MatCopy(appctx->A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
548 ierr = MatScale(A,-1);CHKERRQ(ierr);
549 ierr = MatShift(A,s);CHKERRQ(ierr);
550 ierr = MatCopy(A,B,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
551 appctx->oshift = s;
552 return 0;
553 }
554
555 /*TEST
556
557 test:
558 args: -nox -ts_type ssp -ts_dt 0.0005
559
560 test:
561 suffix: 2
562 args: -nox -ts_type ssp -ts_dt 0.0005 -time_dependent_rhs 1
563
564 test:
565 suffix: 3
566 args: -nox -ts_type rosw -ts_max_steps 3 -ksp_converged_reason
567 filter: sed "s/ATOL/RTOL/g"
568 requires: !single
569
570 test:
571 suffix: 4
572 args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason
573 filter: sed "s/ATOL/RTOL/g"
574
575 test:
576 suffix: 5
577 args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason -time_dependent_rhs
578 filter: sed "s/ATOL/RTOL/g"
579
580 test:
581 requires: !single
582 suffix: pod_guess
583 args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type pod -pc_type none -ksp_converged_reason
584
585 test:
586 requires: !single
587 suffix: pod_guess_Ainner
588 args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type pod -ksp_guess_pod_Ainner -pc_type none -ksp_converged_reason
589
590 test:
591 requires: !single
592 suffix: fischer_guess
593 args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type fischer -pc_type none -ksp_converged_reason
594
595 test:
596 requires: !single
597 suffix: fischer_guess_2
598 args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type fischer -ksp_guess_fischer_model 2,10 -pc_type none -ksp_converged_reason
599
600 test:
601 requires: !single
602 suffix: stringview
603 args: -nox -ts_type rosw -test_string_viewer
604
605 test:
606 requires: !single
607 suffix: stringview_euler
608 args: -nox -ts_type euler -test_string_viewer
609
610 TEST*/
611