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 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);CHKERRMPI(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 */ 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 */ 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 */ 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 */ 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 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 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