1 2 #include <private/tsimpl.h> /*I "petscts.h" I*/ 3 4 /* Logging support */ 5 PetscClassId TS_CLASSID; 6 PetscLogEvent TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval; 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "TSSetTypeFromOptions" 10 /* 11 TSSetTypeFromOptions - Sets the type of ts from user options. 12 13 Collective on TS 14 15 Input Parameter: 16 . ts - The ts 17 18 Level: intermediate 19 20 .keywords: TS, set, options, database, type 21 .seealso: TSSetFromOptions(), TSSetType() 22 */ 23 static PetscErrorCode TSSetTypeFromOptions(TS ts) 24 { 25 PetscBool opt; 26 const char *defaultType; 27 char typeName[256]; 28 PetscErrorCode ierr; 29 30 PetscFunctionBegin; 31 if (((PetscObject)ts)->type_name) { 32 defaultType = ((PetscObject)ts)->type_name; 33 } else { 34 defaultType = TSEULER; 35 } 36 37 if (!TSRegisterAllCalled) {ierr = TSRegisterAll(PETSC_NULL);CHKERRQ(ierr);} 38 ierr = PetscOptionsList("-ts_type", "TS method"," TSSetType", TSList, defaultType, typeName, 256, &opt);CHKERRQ(ierr); 39 if (opt) { 40 ierr = TSSetType(ts, typeName);CHKERRQ(ierr); 41 } else { 42 ierr = TSSetType(ts, defaultType);CHKERRQ(ierr); 43 } 44 PetscFunctionReturn(0); 45 } 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "TSSetFromOptions" 49 /*@ 50 TSSetFromOptions - Sets various TS parameters from user options. 51 52 Collective on TS 53 54 Input Parameter: 55 . ts - the TS context obtained from TSCreate() 56 57 Options Database Keys: 58 + -ts_type <type> - TSEULER, TSBEULER, TSSUNDIALS, TSPSEUDO, TSCN, TSRK, TSTHETA, TSGL, TSSSP 59 . -ts_max_steps maxsteps - maximum number of time-steps to take 60 . -ts_max_time time - maximum time to compute to 61 . -ts_dt dt - initial time step 62 . -ts_monitor - print information at each timestep 63 - -ts_monitor_draw - plot information at each timestep 64 65 Level: beginner 66 67 .keywords: TS, timestep, set, options, database 68 69 .seealso: TSGetType() 70 @*/ 71 PetscErrorCode TSSetFromOptions(TS ts) 72 { 73 PetscReal dt; 74 PetscBool opt,flg; 75 PetscErrorCode ierr; 76 PetscViewer monviewer; 77 char monfilename[PETSC_MAX_PATH_LEN]; 78 SNES snes; 79 80 PetscFunctionBegin; 81 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 82 ierr = PetscOptionsBegin(((PetscObject)ts)->comm, ((PetscObject)ts)->prefix, "Time step options", "TS");CHKERRQ(ierr); 83 /* Handle TS type options */ 84 ierr = TSSetTypeFromOptions(ts);CHKERRQ(ierr); 85 86 /* Handle generic TS options */ 87 ierr = PetscOptionsInt("-ts_max_steps","Maximum number of time steps","TSSetDuration",ts->max_steps,&ts->max_steps,PETSC_NULL);CHKERRQ(ierr); 88 ierr = PetscOptionsReal("-ts_max_time","Time to run to","TSSetDuration",ts->max_time,&ts->max_time,PETSC_NULL);CHKERRQ(ierr); 89 ierr = PetscOptionsReal("-ts_init_time","Initial time","TSSetInitialTime", ts->ptime, &ts->ptime, PETSC_NULL);CHKERRQ(ierr); 90 ierr = PetscOptionsReal("-ts_dt","Initial time step","TSSetInitialTimeStep",ts->initial_time_step,&dt,&opt);CHKERRQ(ierr); 91 if (opt) {ierr = TSSetInitialTimeStep(ts,ts->ptime,dt);CHKERRQ(ierr);} 92 opt = ts->exact_final_time == PETSC_DECIDE ? PETSC_FALSE : (PetscBool)ts->exact_final_time; 93 ierr = PetscOptionsBool("-ts_exact_final_time","Interpolate output to stop exactly at the final time","TSSetExactFinalTime",opt,&opt,&flg);CHKERRQ(ierr); 94 if (flg) {ierr = TSSetExactFinalTime(ts,flg);CHKERRQ(ierr);} 95 ierr = PetscOptionsInt("-ts_max_snes_failures","Maximum number of nonlinear solve failures","",ts->max_snes_failures,&ts->max_snes_failures,PETSC_NULL);CHKERRQ(ierr); 96 ierr = PetscOptionsInt("-ts_max_reject","Maximum number of step rejections","",ts->max_reject,&ts->max_reject,PETSC_NULL);CHKERRQ(ierr); 97 ierr = PetscOptionsBool("-ts_error_if_step_failed","Error if no step succeeds","",ts->errorifstepfailed,&ts->errorifstepfailed,PETSC_NULL);CHKERRQ(ierr); 98 99 /* Monitor options */ 100 ierr = PetscOptionsString("-ts_monitor","Monitor timestep size","TSMonitorDefault","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 101 if (flg) { 102 ierr = PetscViewerASCIIOpen(((PetscObject)ts)->comm,monfilename,&monviewer);CHKERRQ(ierr); 103 ierr = TSMonitorSet(ts,TSMonitorDefault,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);CHKERRQ(ierr); 104 } 105 ierr = PetscOptionsString("-ts_monitor_python","Use Python function","TSMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 106 if (flg) {ierr = PetscPythonMonitorSet((PetscObject)ts,monfilename);CHKERRQ(ierr);} 107 108 opt = PETSC_FALSE; 109 ierr = PetscOptionsBool("-ts_monitor_draw","Monitor timestep size graphically","TSMonitorLG",opt,&opt,PETSC_NULL);CHKERRQ(ierr); 110 if (opt) { 111 ierr = TSMonitorSet(ts,TSMonitorLG,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 112 } 113 opt = PETSC_FALSE; 114 ierr = PetscOptionsBool("-ts_monitor_solution","Monitor solution graphically","TSMonitorSolution",opt,&opt,PETSC_NULL);CHKERRQ(ierr); 115 if (opt) { 116 ierr = TSMonitorSet(ts,TSMonitorSolution,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 117 } 118 119 /* Handle specific TS options */ 120 if (ts->ops->setfromoptions) { 121 ierr = (*ts->ops->setfromoptions)(ts);CHKERRQ(ierr); 122 } 123 124 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 125 ierr = PetscObjectProcessOptionsHandlers((PetscObject)ts);CHKERRQ(ierr); 126 ierr = PetscOptionsEnd();CHKERRQ(ierr); 127 128 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 129 /* Handle subobject options */ 130 if (ts->problem_type == TS_LINEAR) {ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr);} 131 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 132 PetscFunctionReturn(0); 133 } 134 135 #undef __FUNCT__ 136 #undef __FUNCT__ 137 #define __FUNCT__ "TSComputeRHSJacobian" 138 /*@ 139 TSComputeRHSJacobian - Computes the Jacobian matrix that has been 140 set with TSSetRHSJacobian(). 141 142 Collective on TS and Vec 143 144 Input Parameters: 145 + ts - the TS context 146 . t - current timestep 147 - x - input vector 148 149 Output Parameters: 150 + A - Jacobian matrix 151 . B - optional preconditioning matrix 152 - flag - flag indicating matrix structure 153 154 Notes: 155 Most users should not need to explicitly call this routine, as it 156 is used internally within the nonlinear solvers. 157 158 See KSPSetOperators() for important information about setting the 159 flag parameter. 160 161 Level: developer 162 163 .keywords: SNES, compute, Jacobian, matrix 164 165 .seealso: TSSetRHSJacobian(), KSPSetOperators() 166 @*/ 167 PetscErrorCode TSComputeRHSJacobian(TS ts,PetscReal t,Vec X,Mat *A,Mat *B,MatStructure *flg) 168 { 169 PetscErrorCode ierr; 170 PetscInt Xstate; 171 172 PetscFunctionBegin; 173 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 174 PetscValidHeaderSpecific(X,VEC_CLASSID,3); 175 PetscCheckSameComm(ts,1,X,3); 176 ierr = PetscObjectStateQuery((PetscObject)X,&Xstate);CHKERRQ(ierr); 177 if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.X == X && ts->rhsjacobian.Xstate == Xstate))) { 178 *flg = ts->rhsjacobian.mstructure; 179 PetscFunctionReturn(0); 180 } 181 182 if (!ts->userops->rhsjacobian && !ts->userops->ijacobian) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()"); 183 184 if (ts->userops->rhsjacobian) { 185 ierr = PetscLogEventBegin(TS_JacobianEval,ts,X,*A,*B);CHKERRQ(ierr); 186 *flg = DIFFERENT_NONZERO_PATTERN; 187 PetscStackPush("TS user Jacobian function"); 188 ierr = (*ts->userops->rhsjacobian)(ts,t,X,A,B,flg,ts->jacP);CHKERRQ(ierr); 189 PetscStackPop; 190 ierr = PetscLogEventEnd(TS_JacobianEval,ts,X,*A,*B);CHKERRQ(ierr); 191 /* make sure user returned a correct Jacobian and preconditioner */ 192 PetscValidHeaderSpecific(*A,MAT_CLASSID,4); 193 PetscValidHeaderSpecific(*B,MAT_CLASSID,5); 194 } else { 195 ierr = MatZeroEntries(*A);CHKERRQ(ierr); 196 if (*A != *B) {ierr = MatZeroEntries(*B);CHKERRQ(ierr);} 197 *flg = SAME_NONZERO_PATTERN; 198 } 199 ts->rhsjacobian.time = t; 200 ts->rhsjacobian.X = X; 201 ierr = PetscObjectStateQuery((PetscObject)X,&ts->rhsjacobian.Xstate);CHKERRQ(ierr); 202 ts->rhsjacobian.mstructure = *flg; 203 PetscFunctionReturn(0); 204 } 205 206 #undef __FUNCT__ 207 #define __FUNCT__ "TSComputeRHSFunction" 208 /*@ 209 TSComputeRHSFunction - Evaluates the right-hand-side function. 210 211 Collective on TS and Vec 212 213 Input Parameters: 214 + ts - the TS context 215 . t - current time 216 - x - state vector 217 218 Output Parameter: 219 . y - right hand side 220 221 Note: 222 Most users should not need to explicitly call this routine, as it 223 is used internally within the nonlinear solvers. 224 225 Level: developer 226 227 .keywords: TS, compute 228 229 .seealso: TSSetRHSFunction(), TSComputeIFunction() 230 @*/ 231 PetscErrorCode TSComputeRHSFunction(TS ts,PetscReal t,Vec x,Vec y) 232 { 233 PetscErrorCode ierr; 234 235 PetscFunctionBegin; 236 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 237 PetscValidHeaderSpecific(x,VEC_CLASSID,3); 238 PetscValidHeaderSpecific(y,VEC_CLASSID,4); 239 240 if (!ts->userops->rhsfunction && !ts->userops->ifunction) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()"); 241 242 ierr = PetscLogEventBegin(TS_FunctionEval,ts,x,y,0);CHKERRQ(ierr); 243 if (ts->userops->rhsfunction) { 244 PetscStackPush("TS user right-hand-side function"); 245 ierr = (*ts->userops->rhsfunction)(ts,t,x,y,ts->funP);CHKERRQ(ierr); 246 PetscStackPop; 247 } else { 248 ierr = VecZeroEntries(y);CHKERRQ(ierr); 249 } 250 251 ierr = PetscLogEventEnd(TS_FunctionEval,ts,x,y,0);CHKERRQ(ierr); 252 PetscFunctionReturn(0); 253 } 254 255 #undef __FUNCT__ 256 #define __FUNCT__ "TSGetRHSVec_Private" 257 static PetscErrorCode TSGetRHSVec_Private(TS ts,Vec *Frhs) 258 { 259 Vec F; 260 PetscErrorCode ierr; 261 262 PetscFunctionBegin; 263 *Frhs = PETSC_NULL; 264 ierr = TSGetIFunction(ts,&F,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 265 if (!ts->Frhs) { 266 ierr = VecDuplicate(F,&ts->Frhs);CHKERRQ(ierr); 267 } 268 *Frhs = ts->Frhs; 269 PetscFunctionReturn(0); 270 } 271 272 #undef __FUNCT__ 273 #define __FUNCT__ "TSGetRHSMats_Private" 274 static PetscErrorCode TSGetRHSMats_Private(TS ts,Mat *Arhs,Mat *Brhs) 275 { 276 Mat A,B; 277 PetscErrorCode ierr; 278 279 PetscFunctionBegin; 280 ierr = TSGetIJacobian(ts,&A,&B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 281 if (Arhs) { 282 if (!ts->Arhs) { 283 ierr = MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&ts->Arhs);CHKERRQ(ierr); 284 } 285 *Arhs = ts->Arhs; 286 } 287 if (Brhs) { 288 if (!ts->Brhs) { 289 ierr = MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&ts->Brhs);CHKERRQ(ierr); 290 } 291 *Brhs = ts->Brhs; 292 } 293 PetscFunctionReturn(0); 294 } 295 296 #undef __FUNCT__ 297 #define __FUNCT__ "TSComputeIFunction" 298 /*@ 299 TSComputeIFunction - Evaluates the DAE residual written in implicit form F(t,X,Xdot)=0 300 301 Collective on TS and Vec 302 303 Input Parameters: 304 + ts - the TS context 305 . t - current time 306 . X - state vector 307 . Xdot - time derivative of state vector 308 - imex - flag indicates if the method is IMEX so that the RHSFunction should be kept separate 309 310 Output Parameter: 311 . Y - right hand side 312 313 Note: 314 Most users should not need to explicitly call this routine, as it 315 is used internally within the nonlinear solvers. 316 317 If the user did did not write their equations in implicit form, this 318 function recasts them in implicit form. 319 320 Level: developer 321 322 .keywords: TS, compute 323 324 .seealso: TSSetIFunction(), TSComputeRHSFunction() 325 @*/ 326 PetscErrorCode TSComputeIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec Y,PetscBool imex) 327 { 328 PetscErrorCode ierr; 329 330 PetscFunctionBegin; 331 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 332 PetscValidHeaderSpecific(X,VEC_CLASSID,3); 333 PetscValidHeaderSpecific(Xdot,VEC_CLASSID,4); 334 PetscValidHeaderSpecific(Y,VEC_CLASSID,5); 335 336 if (!ts->userops->rhsfunction && !ts->userops->ifunction) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()"); 337 338 ierr = PetscLogEventBegin(TS_FunctionEval,ts,X,Xdot,Y);CHKERRQ(ierr); 339 if (ts->userops->ifunction) { 340 PetscStackPush("TS user implicit function"); 341 ierr = (*ts->userops->ifunction)(ts,t,X,Xdot,Y,ts->funP);CHKERRQ(ierr); 342 PetscStackPop; 343 } 344 if (imex) { 345 if (!ts->userops->ifunction) { 346 ierr = VecCopy(Xdot,Y);CHKERRQ(ierr); 347 } 348 } else if (ts->userops->rhsfunction) { 349 if (ts->userops->ifunction) { 350 Vec Frhs; 351 ierr = TSGetRHSVec_Private(ts,&Frhs);CHKERRQ(ierr); 352 ierr = TSComputeRHSFunction(ts,t,X,Frhs);CHKERRQ(ierr); 353 ierr = VecAXPY(Y,-1,Frhs);CHKERRQ(ierr); 354 } else { 355 ierr = TSComputeRHSFunction(ts,t,X,Y);CHKERRQ(ierr); 356 ierr = VecAYPX(Y,-1,Xdot);CHKERRQ(ierr); 357 } 358 } 359 ierr = PetscLogEventEnd(TS_FunctionEval,ts,X,Xdot,Y);CHKERRQ(ierr); 360 PetscFunctionReturn(0); 361 } 362 363 #undef __FUNCT__ 364 #define __FUNCT__ "TSComputeIJacobian" 365 /*@ 366 TSComputeIJacobian - Evaluates the Jacobian of the DAE 367 368 Collective on TS and Vec 369 370 Input 371 Input Parameters: 372 + ts - the TS context 373 . t - current timestep 374 . X - state vector 375 . Xdot - time derivative of state vector 376 . shift - shift to apply, see note below 377 - imex - flag indicates if the method is IMEX so that the RHSJacobian should be kept separate 378 379 Output Parameters: 380 + A - Jacobian matrix 381 . B - optional preconditioning matrix 382 - flag - flag indicating matrix structure 383 384 Notes: 385 If F(t,X,Xdot)=0 is the DAE, the required Jacobian is 386 387 dF/dX + shift*dF/dXdot 388 389 Most users should not need to explicitly call this routine, as it 390 is used internally within the nonlinear solvers. 391 392 Level: developer 393 394 .keywords: TS, compute, Jacobian, matrix 395 396 .seealso: TSSetIJacobian() 397 @*/ 398 PetscErrorCode TSComputeIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,PetscBool imex) 399 { 400 PetscInt Xstate, Xdotstate; 401 PetscErrorCode ierr; 402 403 PetscFunctionBegin; 404 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 405 PetscValidHeaderSpecific(X,VEC_CLASSID,3); 406 PetscValidHeaderSpecific(Xdot,VEC_CLASSID,4); 407 PetscValidPointer(A,6); 408 PetscValidHeaderSpecific(*A,MAT_CLASSID,6); 409 PetscValidPointer(B,7); 410 PetscValidHeaderSpecific(*B,MAT_CLASSID,7); 411 PetscValidPointer(flg,8); 412 ierr = PetscObjectStateQuery((PetscObject)X,&Xstate);CHKERRQ(ierr); 413 ierr = PetscObjectStateQuery((PetscObject)Xdot,&Xdotstate);CHKERRQ(ierr); 414 if (ts->ijacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->ijacobian.X == X && ts->ijacobian.Xstate == Xstate && ts->ijacobian.Xdot == Xdot && ts->ijacobian.Xdotstate == Xdotstate && ts->ijacobian.imex == imex))) { 415 *flg = ts->ijacobian.mstructure; 416 ierr = MatScale(*A, shift / ts->ijacobian.shift);CHKERRQ(ierr); 417 PetscFunctionReturn(0); 418 } 419 420 if (!ts->userops->rhsjacobian && !ts->userops->ijacobian) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()"); 421 422 *flg = SAME_NONZERO_PATTERN; /* In case we're solving a linear problem in which case it wouldn't get initialized below. */ 423 ierr = PetscLogEventBegin(TS_JacobianEval,ts,X,*A,*B);CHKERRQ(ierr); 424 if (ts->userops->ijacobian) { 425 *flg = DIFFERENT_NONZERO_PATTERN; 426 PetscStackPush("TS user implicit Jacobian"); 427 ierr = (*ts->userops->ijacobian)(ts,t,X,Xdot,shift,A,B,flg,ts->jacP);CHKERRQ(ierr); 428 PetscStackPop; 429 /* make sure user returned a correct Jacobian and preconditioner */ 430 PetscValidHeaderSpecific(*A,MAT_CLASSID,4); 431 PetscValidHeaderSpecific(*B,MAT_CLASSID,5); 432 } 433 if (imex) { 434 if (!ts->userops->ijacobian) { /* system was written as Xdot = F(t,X) */ 435 ierr = MatZeroEntries(*A);CHKERRQ(ierr); 436 ierr = MatShift(*A,shift);CHKERRQ(ierr); 437 if (*A != *B) { 438 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 439 ierr = MatShift(*B,shift);CHKERRQ(ierr); 440 } 441 *flg = SAME_PRECONDITIONER; 442 } 443 } else { 444 if (!ts->userops->ijacobian) { 445 ierr = TSComputeRHSJacobian(ts,t,X,A,B,flg);CHKERRQ(ierr); 446 ierr = MatScale(*A,-1);CHKERRQ(ierr); 447 ierr = MatShift(*A,shift);CHKERRQ(ierr); 448 if (*A != *B) { 449 ierr = MatScale(*B,-1);CHKERRQ(ierr); 450 ierr = MatShift(*B,shift);CHKERRQ(ierr); 451 } 452 } else if (ts->userops->rhsjacobian) { 453 Mat Arhs,Brhs; 454 MatStructure axpy,flg2 = DIFFERENT_NONZERO_PATTERN; 455 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 456 ierr = TSComputeRHSJacobian(ts,t,X,&Arhs,&Brhs,&flg2);CHKERRQ(ierr); 457 axpy = (*flg == flg2) ? SAME_NONZERO_PATTERN : DIFFERENT_NONZERO_PATTERN; 458 ierr = MatAXPY(*A,-1,Arhs,axpy);CHKERRQ(ierr); 459 if (*A != *B) { 460 ierr = MatAXPY(*B,-1,Brhs,axpy);CHKERRQ(ierr); 461 } 462 *flg = PetscMin(*flg,flg2); 463 } 464 } 465 466 ts->ijacobian.time = t; 467 ts->ijacobian.X = X; 468 ts->ijacobian.Xdot = Xdot; 469 ierr = PetscObjectStateQuery((PetscObject)X,&ts->ijacobian.Xstate);CHKERRQ(ierr); 470 ierr = PetscObjectStateQuery((PetscObject)Xdot,&ts->ijacobian.Xdotstate);CHKERRQ(ierr); 471 ts->ijacobian.shift = shift; 472 ts->ijacobian.imex = imex; 473 ts->ijacobian.mstructure = *flg; 474 ierr = PetscLogEventEnd(TS_JacobianEval,ts,X,*A,*B);CHKERRQ(ierr); 475 PetscFunctionReturn(0); 476 } 477 478 #undef __FUNCT__ 479 #define __FUNCT__ "TSSetRHSFunction" 480 /*@C 481 TSSetRHSFunction - Sets the routine for evaluating the function, 482 F(t,u), where U_t = F(t,u). 483 484 Logically Collective on TS 485 486 Input Parameters: 487 + ts - the TS context obtained from TSCreate() 488 . r - vector to put the computed right hand side (or PETSC_NULL to have it created) 489 . f - routine for evaluating the right-hand-side function 490 - ctx - [optional] user-defined context for private data for the 491 function evaluation routine (may be PETSC_NULL) 492 493 Calling sequence of func: 494 $ func (TS ts,PetscReal t,Vec u,Vec F,void *ctx); 495 496 + t - current timestep 497 . u - input vector 498 . F - function vector 499 - ctx - [optional] user-defined function context 500 501 Important: 502 The user MUST call either this routine or TSSetMatrices(). 503 504 Level: beginner 505 506 .keywords: TS, timestep, set, right-hand-side, function 507 508 .seealso: TSSetMatrices() 509 @*/ 510 PetscErrorCode TSSetRHSFunction(TS ts,Vec r,PetscErrorCode (*f)(TS,PetscReal,Vec,Vec,void*),void *ctx) 511 { 512 PetscErrorCode ierr; 513 SNES snes; 514 515 PetscFunctionBegin; 516 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 517 if (r) PetscValidHeaderSpecific(r,VEC_CLASSID,2); 518 if (f) ts->userops->rhsfunction = f; 519 if (ctx) ts->funP = ctx; 520 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 521 ierr = SNESSetFunction(snes,r,SNESTSFormFunction,ts);CHKERRQ(ierr); 522 PetscFunctionReturn(0); 523 } 524 525 #undef __FUNCT__ 526 #define __FUNCT__ "TSSetRHSJacobian" 527 /*@C 528 TSSetRHSJacobian - Sets the function to compute the Jacobian of F, 529 where U_t = F(U,t), as well as the location to store the matrix. 530 Use TSSetMatrices() for linear problems. 531 532 Logically Collective on TS 533 534 Input Parameters: 535 + ts - the TS context obtained from TSCreate() 536 . A - Jacobian matrix 537 . B - preconditioner matrix (usually same as A) 538 . f - the Jacobian evaluation routine 539 - ctx - [optional] user-defined context for private data for the 540 Jacobian evaluation routine (may be PETSC_NULL) 541 542 Calling sequence of func: 543 $ func (TS ts,PetscReal t,Vec u,Mat *A,Mat *B,MatStructure *flag,void *ctx); 544 545 + t - current timestep 546 . u - input vector 547 . A - matrix A, where U_t = A(t)u 548 . B - preconditioner matrix, usually the same as A 549 . flag - flag indicating information about the preconditioner matrix 550 structure (same as flag in KSPSetOperators()) 551 - ctx - [optional] user-defined context for matrix evaluation routine 552 553 Notes: 554 See KSPSetOperators() for important information about setting the flag 555 output parameter in the routine func(). Be sure to read this information! 556 557 The routine func() takes Mat * as the matrix arguments rather than Mat. 558 This allows the matrix evaluation routine to replace A and/or B with a 559 completely new matrix structure (not just different matrix elements) 560 when appropriate, for instance, if the nonzero structure is changing 561 throughout the global iterations. 562 563 Level: beginner 564 565 .keywords: TS, timestep, set, right-hand-side, Jacobian 566 567 .seealso: TSDefaultComputeJacobianColor(), 568 SNESDefaultComputeJacobianColor(), TSSetRHSFunction(), TSSetMatrices() 569 570 @*/ 571 PetscErrorCode TSSetRHSJacobian(TS ts,Mat A,Mat B,TSRHSJacobian f,void *ctx) 572 { 573 PetscErrorCode ierr; 574 SNES snes; 575 576 PetscFunctionBegin; 577 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 578 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,2); 579 if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3); 580 if (A) PetscCheckSameComm(ts,1,A,2); 581 if (B) PetscCheckSameComm(ts,1,B,3); 582 583 if (f) ts->userops->rhsjacobian = f; 584 if (ctx) ts->jacP = ctx; 585 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 586 if (!ts->userops->ijacobian) { 587 ierr = SNESSetJacobian(snes,A,B,SNESTSFormJacobian,ts);CHKERRQ(ierr); 588 } 589 if (A) { 590 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 591 ierr = MatDestroy(&ts->Arhs);CHKERRQ(ierr); 592 ts->Arhs = A; 593 } 594 if (B) { 595 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr); 596 ierr = MatDestroy(&ts->Brhs);CHKERRQ(ierr); 597 ts->Brhs = B; 598 } 599 PetscFunctionReturn(0); 600 } 601 602 603 #undef __FUNCT__ 604 #define __FUNCT__ "TSSetIFunction" 605 /*@C 606 TSSetIFunction - Set the function to compute F(t,U,U_t) where F = 0 is the DAE to be solved. 607 608 Logically Collective on TS 609 610 Input Parameters: 611 + ts - the TS context obtained from TSCreate() 612 . r - vector to hold the residual (or PETSC_NULL to have it created internally) 613 . f - the function evaluation routine 614 - ctx - user-defined context for private data for the function evaluation routine (may be PETSC_NULL) 615 616 Calling sequence of f: 617 $ f(TS ts,PetscReal t,Vec u,Vec u_t,Vec F,ctx); 618 619 + t - time at step/stage being solved 620 . u - state vector 621 . u_t - time derivative of state vector 622 . F - function vector 623 - ctx - [optional] user-defined context for matrix evaluation routine 624 625 Important: 626 The user MUST call either this routine, TSSetRHSFunction(), or TSSetMatrices(). This routine must be used when not solving an ODE. 627 628 Level: beginner 629 630 .keywords: TS, timestep, set, DAE, Jacobian 631 632 .seealso: TSSetMatrices(), TSSetRHSFunction(), TSSetIJacobian() 633 @*/ 634 PetscErrorCode TSSetIFunction(TS ts,Vec res,TSIFunction f,void *ctx) 635 { 636 PetscErrorCode ierr; 637 SNES snes; 638 639 PetscFunctionBegin; 640 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 641 if (res) PetscValidHeaderSpecific(res,VEC_CLASSID,2); 642 if (f) ts->userops->ifunction = f; 643 if (ctx) ts->funP = ctx; 644 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 645 ierr = SNESSetFunction(snes,res,SNESTSFormFunction,ts);CHKERRQ(ierr); 646 PetscFunctionReturn(0); 647 } 648 649 #undef __FUNCT__ 650 #define __FUNCT__ "TSGetIFunction" 651 /*@C 652 TSGetIFunction - Returns the vector where the implicit residual is stored and the function/contex to compute it. 653 654 Not Collective 655 656 Input Parameter: 657 . ts - the TS context 658 659 Output Parameter: 660 + r - vector to hold residual (or PETSC_NULL) 661 . func - the function to compute residual (or PETSC_NULL) 662 - ctx - the function context (or PETSC_NULL) 663 664 Level: advanced 665 666 .keywords: TS, nonlinear, get, function 667 668 .seealso: TSSetIFunction(), SNESGetFunction() 669 @*/ 670 PetscErrorCode TSGetIFunction(TS ts,Vec *r,TSIFunction *func,void **ctx) 671 { 672 PetscErrorCode ierr; 673 SNES snes; 674 675 PetscFunctionBegin; 676 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 677 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 678 ierr = SNESGetFunction(snes,r,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 679 if (func) *func = ts->userops->ifunction; 680 if (ctx) *ctx = ts->funP; 681 PetscFunctionReturn(0); 682 } 683 684 #undef __FUNCT__ 685 #define __FUNCT__ "TSGetRHSFunction" 686 /*@C 687 TSGetRHSFunction - Returns the vector where the right hand side is stored and the function/context to compute it. 688 689 Not Collective 690 691 Input Parameter: 692 . ts - the TS context 693 694 Output Parameter: 695 + r - vector to hold computed right hand side (or PETSC_NULL) 696 . func - the function to compute right hand side (or PETSC_NULL) 697 - ctx - the function context (or PETSC_NULL) 698 699 Level: advanced 700 701 .keywords: TS, nonlinear, get, function 702 703 .seealso: TSSetRhsfunction(), SNESGetFunction() 704 @*/ 705 PetscErrorCode TSGetRHSFunction(TS ts,Vec *r,TSRHSFunction *func,void **ctx) 706 { 707 PetscErrorCode ierr; 708 SNES snes; 709 710 PetscFunctionBegin; 711 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 712 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 713 ierr = SNESGetFunction(snes,r,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 714 if (func) *func = ts->userops->rhsfunction; 715 if (ctx) *ctx = ts->funP; 716 PetscFunctionReturn(0); 717 } 718 719 #undef __FUNCT__ 720 #define __FUNCT__ "TSSetIJacobian" 721 /*@C 722 TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function 723 you provided with TSSetIFunction(). 724 725 Logically Collective on TS 726 727 Input Parameters: 728 + ts - the TS context obtained from TSCreate() 729 . A - Jacobian matrix 730 . B - preconditioning matrix for A (may be same as A) 731 . f - the Jacobian evaluation routine 732 - ctx - user-defined context for private data for the Jacobian evaluation routine (may be PETSC_NULL) 733 734 Calling sequence of f: 735 $ f(TS ts,PetscReal t,Vec U,Vec U_t,PetscReal a,Mat *A,Mat *B,MatStructure *flag,void *ctx); 736 737 + t - time at step/stage being solved 738 . U - state vector 739 . U_t - time derivative of state vector 740 . a - shift 741 . A - Jacobian of G(U) = F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t 742 . B - preconditioning matrix for A, may be same as A 743 . flag - flag indicating information about the preconditioner matrix 744 structure (same as flag in KSPSetOperators()) 745 - ctx - [optional] user-defined context for matrix evaluation routine 746 747 Notes: 748 The matrices A and B are exactly the matrices that are used by SNES for the nonlinear solve. 749 750 The matrix dF/dU + a*dF/dU_t you provide turns out to be 751 the Jacobian of G(U) = F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved. 752 The time integrator internally approximates U_t by W+a*U where the positive "shift" 753 a and vector W depend on the integration method, step size, and past states. For example with 754 the backward Euler method a = 1/dt and W = -a*U(previous timestep) so 755 W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt 756 757 Level: beginner 758 759 .keywords: TS, timestep, DAE, Jacobian 760 761 .seealso: TSSetIFunction(), TSSetRHSJacobian() 762 763 @*/ 764 PetscErrorCode TSSetIJacobian(TS ts,Mat A,Mat B,TSIJacobian f,void *ctx) 765 { 766 PetscErrorCode ierr; 767 SNES snes; 768 769 PetscFunctionBegin; 770 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 771 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,2); 772 if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3); 773 if (A) PetscCheckSameComm(ts,1,A,2); 774 if (B) PetscCheckSameComm(ts,1,B,3); 775 if (f) ts->userops->ijacobian = f; 776 if (ctx) ts->jacP = ctx; 777 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 778 ierr = SNESSetJacobian(snes,A,B,SNESTSFormJacobian,ts);CHKERRQ(ierr); 779 PetscFunctionReturn(0); 780 } 781 782 #undef __FUNCT__ 783 #define __FUNCT__ "TSView" 784 /*@C 785 TSView - Prints the TS data structure. 786 787 Collective on TS 788 789 Input Parameters: 790 + ts - the TS context obtained from TSCreate() 791 - viewer - visualization context 792 793 Options Database Key: 794 . -ts_view - calls TSView() at end of TSStep() 795 796 Notes: 797 The available visualization contexts include 798 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 799 - PETSC_VIEWER_STDOUT_WORLD - synchronized standard 800 output where only the first processor opens 801 the file. All other processors send their 802 data to the first processor to print. 803 804 The user can open an alternative visualization context with 805 PetscViewerASCIIOpen() - output to a specified file. 806 807 Level: beginner 808 809 .keywords: TS, timestep, view 810 811 .seealso: PetscViewerASCIIOpen() 812 @*/ 813 PetscErrorCode TSView(TS ts,PetscViewer viewer) 814 { 815 PetscErrorCode ierr; 816 const TSType type; 817 PetscBool iascii,isstring,isundials; 818 819 PetscFunctionBegin; 820 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 821 if (!viewer) { 822 ierr = PetscViewerASCIIGetStdout(((PetscObject)ts)->comm,&viewer);CHKERRQ(ierr); 823 } 824 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 825 PetscCheckSameComm(ts,1,viewer,2); 826 827 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 828 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 829 if (iascii) { 830 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)ts,viewer,"TS Object");CHKERRQ(ierr); 831 if (ts->ops->view) { 832 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 833 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 834 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 835 } 836 ierr = PetscViewerASCIIPrintf(viewer," maximum steps=%D\n",ts->max_steps);CHKERRQ(ierr); 837 ierr = PetscViewerASCIIPrintf(viewer," maximum time=%G\n",ts->max_time);CHKERRQ(ierr); 838 if (ts->problem_type == TS_NONLINEAR) { 839 ierr = PetscViewerASCIIPrintf(viewer," total number of nonlinear solver iterations=%D\n",ts->nonlinear_its);CHKERRQ(ierr); 840 ierr = PetscViewerASCIIPrintf(viewer," total number of nonlinear solve failures=%D\n",ts->max_snes_failures);CHKERRQ(ierr); 841 } 842 ierr = PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",ts->linear_its);CHKERRQ(ierr); 843 ierr = PetscViewerASCIIPrintf(viewer," total number of rejected steps=%D\n",ts->reject);CHKERRQ(ierr); 844 } else if (isstring) { 845 ierr = TSGetType(ts,&type);CHKERRQ(ierr); 846 ierr = PetscViewerStringSPrintf(viewer," %-7.7s",type);CHKERRQ(ierr); 847 } 848 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 849 ierr = PetscTypeCompare((PetscObject)ts,TSSUNDIALS,&isundials);CHKERRQ(ierr); 850 if (!isundials && ts->snes) {ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr);} 851 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 852 PetscFunctionReturn(0); 853 } 854 855 856 #undef __FUNCT__ 857 #define __FUNCT__ "TSSetApplicationContext" 858 /*@ 859 TSSetApplicationContext - Sets an optional user-defined context for 860 the timesteppers. 861 862 Logically Collective on TS 863 864 Input Parameters: 865 + ts - the TS context obtained from TSCreate() 866 - usrP - optional user context 867 868 Level: intermediate 869 870 .keywords: TS, timestep, set, application, context 871 872 .seealso: TSGetApplicationContext() 873 @*/ 874 PetscErrorCode TSSetApplicationContext(TS ts,void *usrP) 875 { 876 PetscFunctionBegin; 877 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 878 ts->user = usrP; 879 PetscFunctionReturn(0); 880 } 881 882 #undef __FUNCT__ 883 #define __FUNCT__ "TSGetApplicationContext" 884 /*@ 885 TSGetApplicationContext - Gets the user-defined context for the 886 timestepper. 887 888 Not Collective 889 890 Input Parameter: 891 . ts - the TS context obtained from TSCreate() 892 893 Output Parameter: 894 . usrP - user context 895 896 Level: intermediate 897 898 .keywords: TS, timestep, get, application, context 899 900 .seealso: TSSetApplicationContext() 901 @*/ 902 PetscErrorCode TSGetApplicationContext(TS ts,void *usrP) 903 { 904 PetscFunctionBegin; 905 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 906 *(void**)usrP = ts->user; 907 PetscFunctionReturn(0); 908 } 909 910 #undef __FUNCT__ 911 #define __FUNCT__ "TSGetTimeStepNumber" 912 /*@ 913 TSGetTimeStepNumber - Gets the current number of timesteps. 914 915 Not Collective 916 917 Input Parameter: 918 . ts - the TS context obtained from TSCreate() 919 920 Output Parameter: 921 . iter - number steps so far 922 923 Level: intermediate 924 925 .keywords: TS, timestep, get, iteration, number 926 @*/ 927 PetscErrorCode TSGetTimeStepNumber(TS ts,PetscInt* iter) 928 { 929 PetscFunctionBegin; 930 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 931 PetscValidIntPointer(iter,2); 932 *iter = ts->steps; 933 PetscFunctionReturn(0); 934 } 935 936 #undef __FUNCT__ 937 #define __FUNCT__ "TSSetInitialTimeStep" 938 /*@ 939 TSSetInitialTimeStep - Sets the initial timestep to be used, 940 as well as the initial time. 941 942 Logically Collective on TS 943 944 Input Parameters: 945 + ts - the TS context obtained from TSCreate() 946 . initial_time - the initial time 947 - time_step - the size of the timestep 948 949 Level: intermediate 950 951 .seealso: TSSetTimeStep(), TSGetTimeStep() 952 953 .keywords: TS, set, initial, timestep 954 @*/ 955 PetscErrorCode TSSetInitialTimeStep(TS ts,PetscReal initial_time,PetscReal time_step) 956 { 957 PetscErrorCode ierr; 958 959 PetscFunctionBegin; 960 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 961 ierr = TSSetTimeStep(ts,time_step);CHKERRQ(ierr); 962 ts->initial_time_step = time_step; 963 ts->ptime = initial_time; 964 PetscFunctionReturn(0); 965 } 966 967 #undef __FUNCT__ 968 #define __FUNCT__ "TSSetTimeStep" 969 /*@ 970 TSSetTimeStep - Allows one to reset the timestep at any time, 971 useful for simple pseudo-timestepping codes. 972 973 Logically Collective on TS 974 975 Input Parameters: 976 + ts - the TS context obtained from TSCreate() 977 - time_step - the size of the timestep 978 979 Level: intermediate 980 981 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 982 983 .keywords: TS, set, timestep 984 @*/ 985 PetscErrorCode TSSetTimeStep(TS ts,PetscReal time_step) 986 { 987 PetscFunctionBegin; 988 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 989 PetscValidLogicalCollectiveReal(ts,time_step,2); 990 ts->time_step = time_step; 991 ts->next_time_step = time_step; 992 PetscFunctionReturn(0); 993 } 994 995 #undef __FUNCT__ 996 #define __FUNCT__ "TSSetExactFinalTime" 997 /*@ 998 TSSetExactFinalTime - Determines whether to interpolate solution to the 999 exact final time requested by the user or just returns it at the final time 1000 it computed. 1001 1002 Logically Collective on TS 1003 1004 Input Parameter: 1005 + ts - the time-step context 1006 - ft - PETSC_TRUE if interpolates, else PETSC_FALSE 1007 1008 Level: beginner 1009 1010 .seealso: TSSetDuration() 1011 @*/ 1012 PetscErrorCode TSSetExactFinalTime(TS ts,PetscBool flg) 1013 { 1014 1015 PetscFunctionBegin; 1016 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1017 ts->exact_final_time = flg; 1018 PetscFunctionReturn(0); 1019 } 1020 1021 #undef __FUNCT__ 1022 #define __FUNCT__ "TSGetTimeStep" 1023 /*@ 1024 TSGetTimeStep - Gets the current timestep size. 1025 1026 Not Collective 1027 1028 Input Parameter: 1029 . ts - the TS context obtained from TSCreate() 1030 1031 Output Parameter: 1032 . dt - the current timestep size 1033 1034 Level: intermediate 1035 1036 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 1037 1038 .keywords: TS, get, timestep 1039 @*/ 1040 PetscErrorCode TSGetTimeStep(TS ts,PetscReal* dt) 1041 { 1042 PetscFunctionBegin; 1043 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1044 PetscValidDoublePointer(dt,2); 1045 *dt = ts->time_step; 1046 PetscFunctionReturn(0); 1047 } 1048 1049 #undef __FUNCT__ 1050 #define __FUNCT__ "TSGetSolution" 1051 /*@ 1052 TSGetSolution - Returns the solution at the present timestep. It 1053 is valid to call this routine inside the function that you are evaluating 1054 in order to move to the new timestep. This vector not changed until 1055 the solution at the next timestep has been calculated. 1056 1057 Not Collective, but Vec returned is parallel if TS is parallel 1058 1059 Input Parameter: 1060 . ts - the TS context obtained from TSCreate() 1061 1062 Output Parameter: 1063 . v - the vector containing the solution 1064 1065 Level: intermediate 1066 1067 .seealso: TSGetTimeStep() 1068 1069 .keywords: TS, timestep, get, solution 1070 @*/ 1071 PetscErrorCode TSGetSolution(TS ts,Vec *v) 1072 { 1073 PetscFunctionBegin; 1074 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1075 PetscValidPointer(v,2); 1076 *v = ts->vec_sol; 1077 PetscFunctionReturn(0); 1078 } 1079 1080 /* ----- Routines to initialize and destroy a timestepper ---- */ 1081 #undef __FUNCT__ 1082 #define __FUNCT__ "TSSetProblemType" 1083 /*@ 1084 TSSetProblemType - Sets the type of problem to be solved. 1085 1086 Not collective 1087 1088 Input Parameters: 1089 + ts - The TS 1090 - type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms 1091 .vb 1092 U_t = A U 1093 U_t = A(t) U 1094 U_t = F(t,U) 1095 .ve 1096 1097 Level: beginner 1098 1099 .keywords: TS, problem type 1100 .seealso: TSSetUp(), TSProblemType, TS 1101 @*/ 1102 PetscErrorCode TSSetProblemType(TS ts, TSProblemType type) 1103 { 1104 PetscErrorCode ierr; 1105 1106 PetscFunctionBegin; 1107 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1108 ts->problem_type = type; 1109 if (type == TS_LINEAR) { 1110 SNES snes; 1111 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1112 ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 1113 } 1114 PetscFunctionReturn(0); 1115 } 1116 1117 #undef __FUNCT__ 1118 #define __FUNCT__ "TSGetProblemType" 1119 /*@C 1120 TSGetProblemType - Gets the type of problem to be solved. 1121 1122 Not collective 1123 1124 Input Parameter: 1125 . ts - The TS 1126 1127 Output Parameter: 1128 . type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms 1129 .vb 1130 M U_t = A U 1131 M(t) U_t = A(t) U 1132 U_t = F(t,U) 1133 .ve 1134 1135 Level: beginner 1136 1137 .keywords: TS, problem type 1138 .seealso: TSSetUp(), TSProblemType, TS 1139 @*/ 1140 PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type) 1141 { 1142 PetscFunctionBegin; 1143 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1144 PetscValidIntPointer(type,2); 1145 *type = ts->problem_type; 1146 PetscFunctionReturn(0); 1147 } 1148 1149 #undef __FUNCT__ 1150 #define __FUNCT__ "TSSetUp" 1151 /*@ 1152 TSSetUp - Sets up the internal data structures for the later use 1153 of a timestepper. 1154 1155 Collective on TS 1156 1157 Input Parameter: 1158 . ts - the TS context obtained from TSCreate() 1159 1160 Notes: 1161 For basic use of the TS solvers the user need not explicitly call 1162 TSSetUp(), since these actions will automatically occur during 1163 the call to TSStep(). However, if one wishes to control this 1164 phase separately, TSSetUp() should be called after TSCreate() 1165 and optional routines of the form TSSetXXX(), but before TSStep(). 1166 1167 Level: advanced 1168 1169 .keywords: TS, timestep, setup 1170 1171 .seealso: TSCreate(), TSStep(), TSDestroy() 1172 @*/ 1173 PetscErrorCode TSSetUp(TS ts) 1174 { 1175 PetscErrorCode ierr; 1176 1177 PetscFunctionBegin; 1178 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1179 if (ts->setupcalled) PetscFunctionReturn(0); 1180 1181 if (!((PetscObject)ts)->type_name) { 1182 ierr = TSSetType(ts,TSEULER);CHKERRQ(ierr); 1183 } 1184 if (ts->exact_final_time == PETSC_DECIDE) ts->exact_final_time = PETSC_FALSE; 1185 1186 if (!ts->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetSolution() first"); 1187 1188 if (ts->ops->setup) { 1189 ierr = (*ts->ops->setup)(ts);CHKERRQ(ierr); 1190 } 1191 1192 ts->setupcalled = PETSC_TRUE; 1193 PetscFunctionReturn(0); 1194 } 1195 1196 #undef __FUNCT__ 1197 #define __FUNCT__ "TSReset" 1198 /*@ 1199 TSReset - Resets a TS context and removes any allocated Vecs and Mats. 1200 1201 Collective on TS 1202 1203 Input Parameter: 1204 . ts - the TS context obtained from TSCreate() 1205 1206 Level: beginner 1207 1208 .keywords: TS, timestep, reset 1209 1210 .seealso: TSCreate(), TSSetup(), TSDestroy() 1211 @*/ 1212 PetscErrorCode TSReset(TS ts) 1213 { 1214 PetscErrorCode ierr; 1215 1216 PetscFunctionBegin; 1217 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1218 if (ts->ops->reset) { 1219 ierr = (*ts->ops->reset)(ts);CHKERRQ(ierr); 1220 } 1221 if (ts->snes) {ierr = SNESReset(ts->snes);CHKERRQ(ierr);} 1222 ierr = MatDestroy(&ts->Arhs);CHKERRQ(ierr); 1223 ierr = MatDestroy(&ts->Brhs);CHKERRQ(ierr); 1224 ierr = VecDestroy(&ts->Frhs);CHKERRQ(ierr); 1225 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 1226 ierr = VecDestroyVecs(ts->nwork,&ts->work);CHKERRQ(ierr); 1227 ts->setupcalled = PETSC_FALSE; 1228 PetscFunctionReturn(0); 1229 } 1230 1231 #undef __FUNCT__ 1232 #define __FUNCT__ "TSDestroy" 1233 /*@ 1234 TSDestroy - Destroys the timestepper context that was created 1235 with TSCreate(). 1236 1237 Collective on TS 1238 1239 Input Parameter: 1240 . ts - the TS context obtained from TSCreate() 1241 1242 Level: beginner 1243 1244 .keywords: TS, timestepper, destroy 1245 1246 .seealso: TSCreate(), TSSetUp(), TSSolve() 1247 @*/ 1248 PetscErrorCode TSDestroy(TS *ts) 1249 { 1250 PetscErrorCode ierr; 1251 1252 PetscFunctionBegin; 1253 if (!*ts) PetscFunctionReturn(0); 1254 PetscValidHeaderSpecific((*ts),TS_CLASSID,1); 1255 if (--((PetscObject)(*ts))->refct > 0) {*ts = 0; PetscFunctionReturn(0);} 1256 1257 ierr = TSReset((*ts));CHKERRQ(ierr); 1258 1259 /* if memory was published with AMS then destroy it */ 1260 ierr = PetscObjectDepublish((*ts));CHKERRQ(ierr); 1261 if ((*ts)->ops->destroy) {ierr = (*(*ts)->ops->destroy)((*ts));CHKERRQ(ierr);} 1262 1263 ierr = SNESDestroy(&(*ts)->snes);CHKERRQ(ierr); 1264 ierr = DMDestroy(&(*ts)->dm);CHKERRQ(ierr); 1265 ierr = TSMonitorCancel((*ts));CHKERRQ(ierr); 1266 1267 ierr = PetscFree((*ts)->userops); 1268 1269 ierr = PetscHeaderDestroy(ts);CHKERRQ(ierr); 1270 PetscFunctionReturn(0); 1271 } 1272 1273 #undef __FUNCT__ 1274 #define __FUNCT__ "TSGetSNES" 1275 /*@ 1276 TSGetSNES - Returns the SNES (nonlinear solver) associated with 1277 a TS (timestepper) context. Valid only for nonlinear problems. 1278 1279 Not Collective, but SNES is parallel if TS is parallel 1280 1281 Input Parameter: 1282 . ts - the TS context obtained from TSCreate() 1283 1284 Output Parameter: 1285 . snes - the nonlinear solver context 1286 1287 Notes: 1288 The user can then directly manipulate the SNES context to set various 1289 options, etc. Likewise, the user can then extract and manipulate the 1290 KSP, KSP, and PC contexts as well. 1291 1292 TSGetSNES() does not work for integrators that do not use SNES; in 1293 this case TSGetSNES() returns PETSC_NULL in snes. 1294 1295 Level: beginner 1296 1297 .keywords: timestep, get, SNES 1298 @*/ 1299 PetscErrorCode TSGetSNES(TS ts,SNES *snes) 1300 { 1301 PetscErrorCode ierr; 1302 1303 PetscFunctionBegin; 1304 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1305 PetscValidPointer(snes,2); 1306 if (!ts->snes) { 1307 ierr = SNESCreate(((PetscObject)ts)->comm,&ts->snes);CHKERRQ(ierr); 1308 ierr = PetscLogObjectParent(ts,ts->snes);CHKERRQ(ierr); 1309 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);CHKERRQ(ierr); 1310 if (ts->problem_type == TS_LINEAR) { 1311 ierr = SNESSetType(ts->snes,SNESKSPONLY);CHKERRQ(ierr); 1312 } 1313 } 1314 *snes = ts->snes; 1315 PetscFunctionReturn(0); 1316 } 1317 1318 #undef __FUNCT__ 1319 #define __FUNCT__ "TSGetKSP" 1320 /*@ 1321 TSGetKSP - Returns the KSP (linear solver) associated with 1322 a TS (timestepper) context. 1323 1324 Not Collective, but KSP is parallel if TS is parallel 1325 1326 Input Parameter: 1327 . ts - the TS context obtained from TSCreate() 1328 1329 Output Parameter: 1330 . ksp - the nonlinear solver context 1331 1332 Notes: 1333 The user can then directly manipulate the KSP context to set various 1334 options, etc. Likewise, the user can then extract and manipulate the 1335 KSP and PC contexts as well. 1336 1337 TSGetKSP() does not work for integrators that do not use KSP; 1338 in this case TSGetKSP() returns PETSC_NULL in ksp. 1339 1340 Level: beginner 1341 1342 .keywords: timestep, get, KSP 1343 @*/ 1344 PetscErrorCode TSGetKSP(TS ts,KSP *ksp) 1345 { 1346 PetscErrorCode ierr; 1347 SNES snes; 1348 1349 PetscFunctionBegin; 1350 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1351 PetscValidPointer(ksp,2); 1352 if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first"); 1353 if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()"); 1354 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1355 ierr = SNESGetKSP(snes,ksp);CHKERRQ(ierr); 1356 PetscFunctionReturn(0); 1357 } 1358 1359 /* ----------- Routines to set solver parameters ---------- */ 1360 1361 #undef __FUNCT__ 1362 #define __FUNCT__ "TSGetDuration" 1363 /*@ 1364 TSGetDuration - Gets the maximum number of timesteps to use and 1365 maximum time for iteration. 1366 1367 Not Collective 1368 1369 Input Parameters: 1370 + ts - the TS context obtained from TSCreate() 1371 . maxsteps - maximum number of iterations to use, or PETSC_NULL 1372 - maxtime - final time to iterate to, or PETSC_NULL 1373 1374 Level: intermediate 1375 1376 .keywords: TS, timestep, get, maximum, iterations, time 1377 @*/ 1378 PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime) 1379 { 1380 PetscFunctionBegin; 1381 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1382 if (maxsteps) { 1383 PetscValidIntPointer(maxsteps,2); 1384 *maxsteps = ts->max_steps; 1385 } 1386 if (maxtime ) { 1387 PetscValidScalarPointer(maxtime,3); 1388 *maxtime = ts->max_time; 1389 } 1390 PetscFunctionReturn(0); 1391 } 1392 1393 #undef __FUNCT__ 1394 #define __FUNCT__ "TSSetDuration" 1395 /*@ 1396 TSSetDuration - Sets the maximum number of timesteps to use and 1397 maximum time for iteration. 1398 1399 Logically Collective on TS 1400 1401 Input Parameters: 1402 + ts - the TS context obtained from TSCreate() 1403 . maxsteps - maximum number of iterations to use 1404 - maxtime - final time to iterate to 1405 1406 Options Database Keys: 1407 . -ts_max_steps <maxsteps> - Sets maxsteps 1408 . -ts_max_time <maxtime> - Sets maxtime 1409 1410 Notes: 1411 The default maximum number of iterations is 5000. Default time is 5.0 1412 1413 Level: intermediate 1414 1415 .keywords: TS, timestep, set, maximum, iterations 1416 1417 .seealso: TSSetExactFinalTime() 1418 @*/ 1419 PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime) 1420 { 1421 PetscFunctionBegin; 1422 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1423 PetscValidLogicalCollectiveInt(ts,maxsteps,2); 1424 PetscValidLogicalCollectiveReal(ts,maxtime,2); 1425 if (maxsteps >= 0) ts->max_steps = maxsteps; 1426 if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime; 1427 PetscFunctionReturn(0); 1428 } 1429 1430 #undef __FUNCT__ 1431 #define __FUNCT__ "TSSetSolution" 1432 /*@ 1433 TSSetSolution - Sets the initial solution vector 1434 for use by the TS routines. 1435 1436 Logically Collective on TS and Vec 1437 1438 Input Parameters: 1439 + ts - the TS context obtained from TSCreate() 1440 - x - the solution vector 1441 1442 Level: beginner 1443 1444 .keywords: TS, timestep, set, solution, initial conditions 1445 @*/ 1446 PetscErrorCode TSSetSolution(TS ts,Vec x) 1447 { 1448 PetscErrorCode ierr; 1449 1450 PetscFunctionBegin; 1451 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1452 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 1453 ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr); 1454 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 1455 ts->vec_sol = x; 1456 PetscFunctionReturn(0); 1457 } 1458 1459 #undef __FUNCT__ 1460 #define __FUNCT__ "TSSetPreStep" 1461 /*@C 1462 TSSetPreStep - Sets the general-purpose function 1463 called once at the beginning of each time step. 1464 1465 Logically Collective on TS 1466 1467 Input Parameters: 1468 + ts - The TS context obtained from TSCreate() 1469 - func - The function 1470 1471 Calling sequence of func: 1472 . func (TS ts); 1473 1474 Level: intermediate 1475 1476 .keywords: TS, timestep 1477 @*/ 1478 PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS)) 1479 { 1480 PetscFunctionBegin; 1481 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1482 ts->ops->prestep = func; 1483 PetscFunctionReturn(0); 1484 } 1485 1486 #undef __FUNCT__ 1487 #define __FUNCT__ "TSPreStep" 1488 /*@C 1489 TSPreStep - Runs the user-defined pre-step function. 1490 1491 Collective on TS 1492 1493 Input Parameters: 1494 . ts - The TS context obtained from TSCreate() 1495 1496 Notes: 1497 TSPreStep() is typically used within time stepping implementations, 1498 so most users would not generally call this routine themselves. 1499 1500 Level: developer 1501 1502 .keywords: TS, timestep 1503 @*/ 1504 PetscErrorCode TSPreStep(TS ts) 1505 { 1506 PetscErrorCode ierr; 1507 1508 PetscFunctionBegin; 1509 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1510 if (ts->ops->prestep) { 1511 PetscStackPush("TS PreStep function"); 1512 ierr = (*ts->ops->prestep)(ts);CHKERRQ(ierr); 1513 PetscStackPop; 1514 } 1515 PetscFunctionReturn(0); 1516 } 1517 1518 #undef __FUNCT__ 1519 #define __FUNCT__ "TSSetPostStep" 1520 /*@C 1521 TSSetPostStep - Sets the general-purpose function 1522 called once at the end of each time step. 1523 1524 Logically Collective on TS 1525 1526 Input Parameters: 1527 + ts - The TS context obtained from TSCreate() 1528 - func - The function 1529 1530 Calling sequence of func: 1531 . func (TS ts); 1532 1533 Level: intermediate 1534 1535 .keywords: TS, timestep 1536 @*/ 1537 PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS)) 1538 { 1539 PetscFunctionBegin; 1540 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1541 ts->ops->poststep = func; 1542 PetscFunctionReturn(0); 1543 } 1544 1545 #undef __FUNCT__ 1546 #define __FUNCT__ "TSPostStep" 1547 /*@C 1548 TSPostStep - Runs the user-defined post-step function. 1549 1550 Collective on TS 1551 1552 Input Parameters: 1553 . ts - The TS context obtained from TSCreate() 1554 1555 Notes: 1556 TSPostStep() is typically used within time stepping implementations, 1557 so most users would not generally call this routine themselves. 1558 1559 Level: developer 1560 1561 .keywords: TS, timestep 1562 @*/ 1563 PetscErrorCode TSPostStep(TS ts) 1564 { 1565 PetscErrorCode ierr; 1566 1567 PetscFunctionBegin; 1568 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1569 if (ts->ops->poststep) { 1570 PetscStackPush("TS PostStep function"); 1571 ierr = (*ts->ops->poststep)(ts);CHKERRQ(ierr); 1572 PetscStackPop; 1573 } 1574 PetscFunctionReturn(0); 1575 } 1576 1577 /* ------------ Routines to set performance monitoring options ----------- */ 1578 1579 #undef __FUNCT__ 1580 #define __FUNCT__ "TSMonitorSet" 1581 /*@C 1582 TSMonitorSet - Sets an ADDITIONAL function that is to be used at every 1583 timestep to display the iteration's progress. 1584 1585 Logically Collective on TS 1586 1587 Input Parameters: 1588 + ts - the TS context obtained from TSCreate() 1589 . func - monitoring routine 1590 . mctx - [optional] user-defined context for private data for the 1591 monitor routine (use PETSC_NULL if no context is desired) 1592 - monitordestroy - [optional] routine that frees monitor context 1593 (may be PETSC_NULL) 1594 1595 Calling sequence of func: 1596 $ int func(TS ts,PetscInt steps,PetscReal time,Vec x,void *mctx) 1597 1598 + ts - the TS context 1599 . steps - iteration number 1600 . time - current time 1601 . x - current iterate 1602 - mctx - [optional] monitoring context 1603 1604 Notes: 1605 This routine adds an additional monitor to the list of monitors that 1606 already has been loaded. 1607 1608 Fortran notes: Only a single monitor function can be set for each TS object 1609 1610 Level: intermediate 1611 1612 .keywords: TS, timestep, set, monitor 1613 1614 .seealso: TSMonitorDefault(), TSMonitorCancel() 1615 @*/ 1616 PetscErrorCode TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**)) 1617 { 1618 PetscFunctionBegin; 1619 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1620 if (ts->numbermonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 1621 ts->monitor[ts->numbermonitors] = monitor; 1622 ts->mdestroy[ts->numbermonitors] = mdestroy; 1623 ts->monitorcontext[ts->numbermonitors++] = (void*)mctx; 1624 PetscFunctionReturn(0); 1625 } 1626 1627 #undef __FUNCT__ 1628 #define __FUNCT__ "TSMonitorCancel" 1629 /*@C 1630 TSMonitorCancel - Clears all the monitors that have been set on a time-step object. 1631 1632 Logically Collective on TS 1633 1634 Input Parameters: 1635 . ts - the TS context obtained from TSCreate() 1636 1637 Notes: 1638 There is no way to remove a single, specific monitor. 1639 1640 Level: intermediate 1641 1642 .keywords: TS, timestep, set, monitor 1643 1644 .seealso: TSMonitorDefault(), TSMonitorSet() 1645 @*/ 1646 PetscErrorCode TSMonitorCancel(TS ts) 1647 { 1648 PetscErrorCode ierr; 1649 PetscInt i; 1650 1651 PetscFunctionBegin; 1652 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1653 for (i=0; i<ts->numbermonitors; i++) { 1654 if (ts->mdestroy[i]) { 1655 ierr = (*ts->mdestroy[i])(&ts->monitorcontext[i]);CHKERRQ(ierr); 1656 } 1657 } 1658 ts->numbermonitors = 0; 1659 PetscFunctionReturn(0); 1660 } 1661 1662 #undef __FUNCT__ 1663 #define __FUNCT__ "TSMonitorDefault" 1664 /*@ 1665 TSMonitorDefault - Sets the Default monitor 1666 1667 Level: intermediate 1668 1669 .keywords: TS, set, monitor 1670 1671 .seealso: TSMonitorDefault(), TSMonitorSet() 1672 @*/ 1673 PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy) 1674 { 1675 PetscErrorCode ierr; 1676 PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(((PetscObject)ts)->comm); 1677 1678 PetscFunctionBegin; 1679 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 1680 ierr = PetscViewerASCIIPrintf(viewer,"%D TS dt %G time %G\n",step,ts->time_step,ptime);CHKERRQ(ierr); 1681 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 1682 PetscFunctionReturn(0); 1683 } 1684 1685 #undef __FUNCT__ 1686 #define __FUNCT__ "TSSetRetainStages" 1687 /*@ 1688 TSSetRetainStages - Request that all stages in the upcoming step be stored so that interpolation will be available. 1689 1690 Logically Collective on TS 1691 1692 Input Argument: 1693 . ts - time stepping context 1694 1695 Output Argument: 1696 . flg - PETSC_TRUE or PETSC_FALSE 1697 1698 Level: intermediate 1699 1700 .keywords: TS, set 1701 1702 .seealso: TSInterpolate(), TSSetPostStep() 1703 @*/ 1704 PetscErrorCode TSSetRetainStages(TS ts,PetscBool flg) 1705 { 1706 1707 PetscFunctionBegin; 1708 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1709 ts->retain_stages = flg; 1710 PetscFunctionReturn(0); 1711 } 1712 1713 #undef __FUNCT__ 1714 #define __FUNCT__ "TSInterpolate" 1715 /*@ 1716 TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval 1717 1718 Collective on TS 1719 1720 Input Argument: 1721 + ts - time stepping context 1722 - t - time to interpolate to 1723 1724 Output Argument: 1725 . X - state at given time 1726 1727 Notes: 1728 The user should call TSSetRetainStages() before taking a step in which interpolation will be requested. 1729 1730 Level: intermediate 1731 1732 Developer Notes: 1733 TSInterpolate() and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints. 1734 1735 .keywords: TS, set 1736 1737 .seealso: TSSetRetainStages(), TSSetPostStep() 1738 @*/ 1739 PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec X) 1740 { 1741 PetscErrorCode ierr; 1742 1743 PetscFunctionBegin; 1744 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1745 if (t < ts->ptime - ts->time_step || ts->ptime < t) SETERRQ3(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Requested time %G not in last time steps [%G,%G]",t,ts->ptime-ts->time_step,ts->ptime); 1746 if (!ts->ops->interpolate) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name); 1747 ierr = (*ts->ops->interpolate)(ts,t,X);CHKERRQ(ierr); 1748 PetscFunctionReturn(0); 1749 } 1750 1751 #undef __FUNCT__ 1752 #define __FUNCT__ "TSStep" 1753 /*@ 1754 TSStep - Steps the requested number of timesteps. 1755 1756 Collective on TS 1757 1758 Input Parameter: 1759 . ts - the TS context obtained from TSCreate() 1760 1761 Output Parameters: 1762 + steps - number of iterations until termination 1763 - ptime - time until termination 1764 1765 Level: beginner 1766 1767 .keywords: TS, timestep, solve 1768 1769 .seealso: TSCreate(), TSSetUp(), TSDestroy() 1770 @*/ 1771 PetscErrorCode TSStep(TS ts) 1772 { 1773 PetscErrorCode ierr; 1774 1775 PetscFunctionBegin; 1776 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1777 1778 ierr = TSSetUp(ts);CHKERRQ(ierr); 1779 1780 ierr = PetscLogEventBegin(TS_Step, ts, 0, 0, 0);CHKERRQ(ierr); 1781 ierr = (*ts->ops->step)(ts);CHKERRQ(ierr); 1782 ierr = PetscLogEventEnd(TS_Step, ts, 0, 0, 0);CHKERRQ(ierr); 1783 PetscFunctionReturn(0); 1784 } 1785 1786 #undef __FUNCT__ 1787 #define __FUNCT__ "TSSolve" 1788 /*@ 1789 TSSolve - Steps the requested number of timesteps. 1790 1791 Collective on TS 1792 1793 Input Parameter: 1794 + ts - the TS context obtained from TSCreate() 1795 - x - the solution vector, or PETSC_NULL if it was set with TSSetSolution() 1796 1797 Output Parameter: 1798 . ftime - time of the state vector x upon completion 1799 1800 Level: beginner 1801 1802 Notes: 1803 The final time returned by this function may be different from the time of the internally 1804 held state accessible by TSGetSolution() and TSGetTime() because the method may have 1805 stepped over the final time. 1806 1807 .keywords: TS, timestep, solve 1808 1809 .seealso: TSCreate(), TSSetSolution(), TSStep() 1810 @*/ 1811 PetscErrorCode TSSolve(TS ts,Vec x,PetscReal *ftime) 1812 { 1813 PetscInt i; 1814 PetscBool flg; 1815 char filename[PETSC_MAX_PATH_LEN]; 1816 PetscViewer viewer; 1817 PetscErrorCode ierr; 1818 1819 PetscFunctionBegin; 1820 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1821 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 1822 if (ts->exact_final_time) { /* Need ts->vec_sol to be distinct so it is not overwritten when we interpolate at the end */ 1823 if (!ts->vec_sol || x == ts->vec_sol) { 1824 Vec y; 1825 ierr = VecDuplicate(x,&y);CHKERRQ(ierr); 1826 ierr = VecCopy(x,y);CHKERRQ(ierr); 1827 ierr = TSSetSolution(ts,y);CHKERRQ(ierr); 1828 ierr = VecDestroy(&y);CHKERRQ(ierr); /* grant ownership */ 1829 } else { 1830 ierr = VecCopy(x,ts->vec_sol);CHKERRQ(ierr); 1831 } 1832 } else { 1833 ierr = TSSetSolution(ts,x);CHKERRQ(ierr); 1834 } 1835 ierr = TSSetUp(ts);CHKERRQ(ierr); 1836 /* reset time step and iteration counters */ 1837 ts->steps = 0; 1838 ts->linear_its = 0; 1839 ts->nonlinear_its = 0; 1840 ts->reason = TS_CONVERGED_ITERATING; 1841 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 1842 1843 if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */ 1844 ierr = (*ts->ops->solve)(ts);CHKERRQ(ierr); 1845 ierr = VecCopy(ts->vec_sol,x);CHKERRQ(ierr); 1846 if (*ftime) *ftime = ts->ptime; 1847 } else { 1848 i = 0; 1849 if (i >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 1850 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 1851 /* steps the requested number of timesteps. */ 1852 while (!ts->reason) { 1853 ierr = TSPreStep(ts);CHKERRQ(ierr); 1854 ierr = TSStep(ts);CHKERRQ(ierr); 1855 if (ts->reason < 0) { 1856 if (ts->errorifstepfailed) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_NOT_CONVERGED,"TSStep has failed"); 1857 } else if (++i >= ts->max_steps) { 1858 ts->reason = TS_CONVERGED_ITS; 1859 } else if (ts->ptime >= ts->max_time) { 1860 ts->reason = TS_CONVERGED_TIME; 1861 } 1862 ierr = TSPostStep(ts);CHKERRQ(ierr); 1863 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 1864 } 1865 if (ts->exact_final_time && ts->ptime >= ts->max_time) { 1866 ierr = TSInterpolate(ts,ts->max_time,x);CHKERRQ(ierr); 1867 if (ftime) *ftime = ts->max_time; 1868 } else { 1869 if (x != ts->vec_sol) { 1870 ierr = VecCopy(ts->vec_sol,x);CHKERRQ(ierr); 1871 } 1872 if (ftime) *ftime = ts->ptime; 1873 } 1874 } 1875 ierr = PetscOptionsGetString(((PetscObject)ts)->prefix,"-ts_view",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 1876 if (flg && !PetscPreLoadingOn) { 1877 ierr = PetscViewerASCIIOpen(((PetscObject)ts)->comm,filename,&viewer);CHKERRQ(ierr); 1878 ierr = TSView(ts,viewer);CHKERRQ(ierr); 1879 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 1880 } 1881 PetscFunctionReturn(0); 1882 } 1883 1884 #undef __FUNCT__ 1885 #define __FUNCT__ "TSMonitor" 1886 /* 1887 Runs the user provided monitor routines, if they exists. 1888 */ 1889 PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec x) 1890 { 1891 PetscErrorCode ierr; 1892 PetscInt i,n = ts->numbermonitors; 1893 1894 PetscFunctionBegin; 1895 for (i=0; i<n; i++) { 1896 ierr = (*ts->monitor[i])(ts,step,ptime,x,ts->monitorcontext[i]);CHKERRQ(ierr); 1897 } 1898 PetscFunctionReturn(0); 1899 } 1900 1901 /* ------------------------------------------------------------------------*/ 1902 1903 #undef __FUNCT__ 1904 #define __FUNCT__ "TSMonitorLGCreate" 1905 /*@C 1906 TSMonitorLGCreate - Creates a line graph context for use with 1907 TS to monitor convergence of preconditioned residual norms. 1908 1909 Collective on TS 1910 1911 Input Parameters: 1912 + host - the X display to open, or null for the local machine 1913 . label - the title to put in the title bar 1914 . x, y - the screen coordinates of the upper left coordinate of the window 1915 - m, n - the screen width and height in pixels 1916 1917 Output Parameter: 1918 . draw - the drawing context 1919 1920 Options Database Key: 1921 . -ts_monitor_draw - automatically sets line graph monitor 1922 1923 Notes: 1924 Use TSMonitorLGDestroy() to destroy this line graph, not PetscDrawLGDestroy(). 1925 1926 Level: intermediate 1927 1928 .keywords: TS, monitor, line graph, residual, seealso 1929 1930 .seealso: TSMonitorLGDestroy(), TSMonitorSet() 1931 1932 @*/ 1933 PetscErrorCode TSMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 1934 { 1935 PetscDraw win; 1936 PetscErrorCode ierr; 1937 1938 PetscFunctionBegin; 1939 ierr = PetscDrawCreate(PETSC_COMM_SELF,host,label,x,y,m,n,&win);CHKERRQ(ierr); 1940 ierr = PetscDrawSetType(win,PETSC_DRAW_X);CHKERRQ(ierr); 1941 ierr = PetscDrawLGCreate(win,1,draw);CHKERRQ(ierr); 1942 ierr = PetscDrawLGIndicateDataPoints(*draw);CHKERRQ(ierr); 1943 1944 ierr = PetscLogObjectParent(*draw,win);CHKERRQ(ierr); 1945 PetscFunctionReturn(0); 1946 } 1947 1948 #undef __FUNCT__ 1949 #define __FUNCT__ "TSMonitorLG" 1950 PetscErrorCode TSMonitorLG(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 1951 { 1952 PetscDrawLG lg = (PetscDrawLG) monctx; 1953 PetscReal x,y = ptime; 1954 PetscErrorCode ierr; 1955 1956 PetscFunctionBegin; 1957 if (!monctx) { 1958 MPI_Comm comm; 1959 PetscViewer viewer; 1960 1961 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 1962 viewer = PETSC_VIEWER_DRAW_(comm); 1963 ierr = PetscViewerDrawGetDrawLG(viewer,0,&lg);CHKERRQ(ierr); 1964 } 1965 1966 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 1967 x = (PetscReal)n; 1968 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 1969 if (n < 20 || (n % 5)) { 1970 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 1971 } 1972 PetscFunctionReturn(0); 1973 } 1974 1975 #undef __FUNCT__ 1976 #define __FUNCT__ "TSMonitorLGDestroy" 1977 /*@C 1978 TSMonitorLGDestroy - Destroys a line graph context that was created 1979 with TSMonitorLGCreate(). 1980 1981 Collective on PetscDrawLG 1982 1983 Input Parameter: 1984 . draw - the drawing context 1985 1986 Level: intermediate 1987 1988 .keywords: TS, monitor, line graph, destroy 1989 1990 .seealso: TSMonitorLGCreate(), TSMonitorSet(), TSMonitorLG(); 1991 @*/ 1992 PetscErrorCode TSMonitorLGDestroy(PetscDrawLG *drawlg) 1993 { 1994 PetscDraw draw; 1995 PetscErrorCode ierr; 1996 1997 PetscFunctionBegin; 1998 ierr = PetscDrawLGGetDraw(*drawlg,&draw);CHKERRQ(ierr); 1999 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 2000 ierr = PetscDrawLGDestroy(drawlg);CHKERRQ(ierr); 2001 PetscFunctionReturn(0); 2002 } 2003 2004 #undef __FUNCT__ 2005 #define __FUNCT__ "TSGetTime" 2006 /*@ 2007 TSGetTime - Gets the current time. 2008 2009 Not Collective 2010 2011 Input Parameter: 2012 . ts - the TS context obtained from TSCreate() 2013 2014 Output Parameter: 2015 . t - the current time 2016 2017 Level: beginner 2018 2019 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 2020 2021 .keywords: TS, get, time 2022 @*/ 2023 PetscErrorCode TSGetTime(TS ts,PetscReal* t) 2024 { 2025 PetscFunctionBegin; 2026 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2027 PetscValidDoublePointer(t,2); 2028 *t = ts->ptime; 2029 PetscFunctionReturn(0); 2030 } 2031 2032 #undef __FUNCT__ 2033 #define __FUNCT__ "TSSetTime" 2034 /*@ 2035 TSSetTime - Allows one to reset the time. 2036 2037 Logically Collective on TS 2038 2039 Input Parameters: 2040 + ts - the TS context obtained from TSCreate() 2041 - time - the time 2042 2043 Level: intermediate 2044 2045 .seealso: TSGetTime(), TSSetDuration() 2046 2047 .keywords: TS, set, time 2048 @*/ 2049 PetscErrorCode TSSetTime(TS ts, PetscReal t) 2050 { 2051 PetscFunctionBegin; 2052 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2053 PetscValidLogicalCollectiveReal(ts,t,2); 2054 ts->ptime = t; 2055 PetscFunctionReturn(0); 2056 } 2057 2058 #undef __FUNCT__ 2059 #define __FUNCT__ "TSSetOptionsPrefix" 2060 /*@C 2061 TSSetOptionsPrefix - Sets the prefix used for searching for all 2062 TS options in the database. 2063 2064 Logically Collective on TS 2065 2066 Input Parameter: 2067 + ts - The TS context 2068 - prefix - The prefix to prepend to all option names 2069 2070 Notes: 2071 A hyphen (-) must NOT be given at the beginning of the prefix name. 2072 The first character of all runtime options is AUTOMATICALLY the 2073 hyphen. 2074 2075 Level: advanced 2076 2077 .keywords: TS, set, options, prefix, database 2078 2079 .seealso: TSSetFromOptions() 2080 2081 @*/ 2082 PetscErrorCode TSSetOptionsPrefix(TS ts,const char prefix[]) 2083 { 2084 PetscErrorCode ierr; 2085 SNES snes; 2086 2087 PetscFunctionBegin; 2088 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2089 ierr = PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2090 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2091 ierr = SNESSetOptionsPrefix(snes,prefix);CHKERRQ(ierr); 2092 PetscFunctionReturn(0); 2093 } 2094 2095 2096 #undef __FUNCT__ 2097 #define __FUNCT__ "TSAppendOptionsPrefix" 2098 /*@C 2099 TSAppendOptionsPrefix - Appends to the prefix used for searching for all 2100 TS options in the database. 2101 2102 Logically Collective on TS 2103 2104 Input Parameter: 2105 + ts - The TS context 2106 - prefix - The prefix to prepend to all option names 2107 2108 Notes: 2109 A hyphen (-) must NOT be given at the beginning of the prefix name. 2110 The first character of all runtime options is AUTOMATICALLY the 2111 hyphen. 2112 2113 Level: advanced 2114 2115 .keywords: TS, append, options, prefix, database 2116 2117 .seealso: TSGetOptionsPrefix() 2118 2119 @*/ 2120 PetscErrorCode TSAppendOptionsPrefix(TS ts,const char prefix[]) 2121 { 2122 PetscErrorCode ierr; 2123 SNES snes; 2124 2125 PetscFunctionBegin; 2126 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2127 ierr = PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2128 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2129 ierr = SNESAppendOptionsPrefix(snes,prefix);CHKERRQ(ierr); 2130 PetscFunctionReturn(0); 2131 } 2132 2133 #undef __FUNCT__ 2134 #define __FUNCT__ "TSGetOptionsPrefix" 2135 /*@C 2136 TSGetOptionsPrefix - Sets the prefix used for searching for all 2137 TS options in the database. 2138 2139 Not Collective 2140 2141 Input Parameter: 2142 . ts - The TS context 2143 2144 Output Parameter: 2145 . prefix - A pointer to the prefix string used 2146 2147 Notes: On the fortran side, the user should pass in a string 'prifix' of 2148 sufficient length to hold the prefix. 2149 2150 Level: intermediate 2151 2152 .keywords: TS, get, options, prefix, database 2153 2154 .seealso: TSAppendOptionsPrefix() 2155 @*/ 2156 PetscErrorCode TSGetOptionsPrefix(TS ts,const char *prefix[]) 2157 { 2158 PetscErrorCode ierr; 2159 2160 PetscFunctionBegin; 2161 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2162 PetscValidPointer(prefix,2); 2163 ierr = PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2164 PetscFunctionReturn(0); 2165 } 2166 2167 #undef __FUNCT__ 2168 #define __FUNCT__ "TSGetRHSJacobian" 2169 /*@C 2170 TSGetRHSJacobian - Returns the Jacobian J at the present timestep. 2171 2172 Not Collective, but parallel objects are returned if TS is parallel 2173 2174 Input Parameter: 2175 . ts - The TS context obtained from TSCreate() 2176 2177 Output Parameters: 2178 + J - The Jacobian J of F, where U_t = F(U,t) 2179 . M - The preconditioner matrix, usually the same as J 2180 . func - Function to compute the Jacobian of the RHS 2181 - ctx - User-defined context for Jacobian evaluation routine 2182 2183 Notes: You can pass in PETSC_NULL for any return argument you do not need. 2184 2185 Level: intermediate 2186 2187 .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 2188 2189 .keywords: TS, timestep, get, matrix, Jacobian 2190 @*/ 2191 PetscErrorCode TSGetRHSJacobian(TS ts,Mat *J,Mat *M,TSRHSJacobian *func,void **ctx) 2192 { 2193 PetscErrorCode ierr; 2194 SNES snes; 2195 2196 PetscFunctionBegin; 2197 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2198 ierr = SNESGetJacobian(snes,J,M,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2199 if (func) *func = ts->userops->rhsjacobian; 2200 if (ctx) *ctx = ts->jacP; 2201 PetscFunctionReturn(0); 2202 } 2203 2204 #undef __FUNCT__ 2205 #define __FUNCT__ "TSGetIJacobian" 2206 /*@C 2207 TSGetIJacobian - Returns the implicit Jacobian at the present timestep. 2208 2209 Not Collective, but parallel objects are returned if TS is parallel 2210 2211 Input Parameter: 2212 . ts - The TS context obtained from TSCreate() 2213 2214 Output Parameters: 2215 + A - The Jacobian of F(t,U,U_t) 2216 . B - The preconditioner matrix, often the same as A 2217 . f - The function to compute the matrices 2218 - ctx - User-defined context for Jacobian evaluation routine 2219 2220 Notes: You can pass in PETSC_NULL for any return argument you do not need. 2221 2222 Level: advanced 2223 2224 .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 2225 2226 .keywords: TS, timestep, get, matrix, Jacobian 2227 @*/ 2228 PetscErrorCode TSGetIJacobian(TS ts,Mat *A,Mat *B,TSIJacobian *f,void **ctx) 2229 { 2230 PetscErrorCode ierr; 2231 SNES snes; 2232 2233 PetscFunctionBegin; 2234 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2235 ierr = SNESGetJacobian(snes,A,B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2236 if (f) *f = ts->userops->ijacobian; 2237 if (ctx) *ctx = ts->jacP; 2238 PetscFunctionReturn(0); 2239 } 2240 2241 #undef __FUNCT__ 2242 #define __FUNCT__ "TSMonitorSolution" 2243 /*@C 2244 TSMonitorSolution - Monitors progress of the TS solvers by calling 2245 VecView() for the solution at each timestep 2246 2247 Collective on TS 2248 2249 Input Parameters: 2250 + ts - the TS context 2251 . step - current time-step 2252 . ptime - current time 2253 - dummy - either a viewer or PETSC_NULL 2254 2255 Level: intermediate 2256 2257 .keywords: TS, vector, monitor, view 2258 2259 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 2260 @*/ 2261 PetscErrorCode TSMonitorSolution(TS ts,PetscInt step,PetscReal ptime,Vec x,void *dummy) 2262 { 2263 PetscErrorCode ierr; 2264 PetscViewer viewer = (PetscViewer) dummy; 2265 2266 PetscFunctionBegin; 2267 if (!dummy) { 2268 viewer = PETSC_VIEWER_DRAW_(((PetscObject)ts)->comm); 2269 } 2270 ierr = VecView(x,viewer);CHKERRQ(ierr); 2271 PetscFunctionReturn(0); 2272 } 2273 2274 2275 #undef __FUNCT__ 2276 #define __FUNCT__ "TSSetDM" 2277 /*@ 2278 TSSetDM - Sets the DM that may be used by some preconditioners 2279 2280 Logically Collective on TS and DM 2281 2282 Input Parameters: 2283 + ts - the preconditioner context 2284 - dm - the dm 2285 2286 Level: intermediate 2287 2288 2289 .seealso: TSGetDM(), SNESSetDM(), SNESGetDM() 2290 @*/ 2291 PetscErrorCode TSSetDM(TS ts,DM dm) 2292 { 2293 PetscErrorCode ierr; 2294 SNES snes; 2295 2296 PetscFunctionBegin; 2297 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2298 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 2299 ierr = DMDestroy(&ts->dm);CHKERRQ(ierr); 2300 ts->dm = dm; 2301 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2302 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 2303 PetscFunctionReturn(0); 2304 } 2305 2306 #undef __FUNCT__ 2307 #define __FUNCT__ "TSGetDM" 2308 /*@ 2309 TSGetDM - Gets the DM that may be used by some preconditioners 2310 2311 Not Collective 2312 2313 Input Parameter: 2314 . ts - the preconditioner context 2315 2316 Output Parameter: 2317 . dm - the dm 2318 2319 Level: intermediate 2320 2321 2322 .seealso: TSSetDM(), SNESSetDM(), SNESGetDM() 2323 @*/ 2324 PetscErrorCode TSGetDM(TS ts,DM *dm) 2325 { 2326 PetscFunctionBegin; 2327 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2328 *dm = ts->dm; 2329 PetscFunctionReturn(0); 2330 } 2331 2332 #undef __FUNCT__ 2333 #define __FUNCT__ "SNESTSFormFunction" 2334 /*@ 2335 SNESTSFormFunction - Function to evaluate nonlinear residual 2336 2337 Logically Collective on SNES 2338 2339 Input Parameter: 2340 + snes - nonlinear solver 2341 . X - the current state at which to evaluate the residual 2342 - ctx - user context, must be a TS 2343 2344 Output Parameter: 2345 . F - the nonlinear residual 2346 2347 Notes: 2348 This function is not normally called by users and is automatically registered with the SNES used by TS. 2349 It is most frequently passed to MatFDColoringSetFunction(). 2350 2351 Level: advanced 2352 2353 .seealso: SNESSetFunction(), MatFDColoringSetFunction() 2354 @*/ 2355 PetscErrorCode SNESTSFormFunction(SNES snes,Vec X,Vec F,void *ctx) 2356 { 2357 TS ts = (TS)ctx; 2358 PetscErrorCode ierr; 2359 2360 PetscFunctionBegin; 2361 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2362 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2363 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 2364 PetscValidHeaderSpecific(ts,TS_CLASSID,4); 2365 ierr = (ts->ops->snesfunction)(snes,X,F,ts);CHKERRQ(ierr); 2366 PetscFunctionReturn(0); 2367 } 2368 2369 #undef __FUNCT__ 2370 #define __FUNCT__ "SNESTSFormJacobian" 2371 /*@ 2372 SNESTSFormJacobian - Function to evaluate the Jacobian 2373 2374 Collective on SNES 2375 2376 Input Parameter: 2377 + snes - nonlinear solver 2378 . X - the current state at which to evaluate the residual 2379 - ctx - user context, must be a TS 2380 2381 Output Parameter: 2382 + A - the Jacobian 2383 . B - the preconditioning matrix (may be the same as A) 2384 - flag - indicates any structure change in the matrix 2385 2386 Notes: 2387 This function is not normally called by users and is automatically registered with the SNES used by TS. 2388 2389 Level: developer 2390 2391 .seealso: SNESSetJacobian() 2392 @*/ 2393 PetscErrorCode SNESTSFormJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flag,void *ctx) 2394 { 2395 TS ts = (TS)ctx; 2396 PetscErrorCode ierr; 2397 2398 PetscFunctionBegin; 2399 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2400 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2401 PetscValidPointer(A,3); 2402 PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 2403 PetscValidPointer(B,4); 2404 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); 2405 PetscValidPointer(flag,5); 2406 PetscValidHeaderSpecific(ts,TS_CLASSID,6); 2407 ierr = (ts->ops->snesjacobian)(snes,X,A,B,flag,ts);CHKERRQ(ierr); 2408 PetscFunctionReturn(0); 2409 } 2410 2411 #undef __FUNCT__ 2412 #define __FUNCT__ "TSComputeRHSFunctionLinear" 2413 /*@C 2414 TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems only 2415 2416 Collective on TS 2417 2418 Input Arguments: 2419 + ts - time stepping context 2420 . t - time at which to evaluate 2421 . X - state at which to evaluate 2422 - ctx - context 2423 2424 Output Arguments: 2425 . F - right hand side 2426 2427 Level: intermediate 2428 2429 Notes: 2430 This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems. 2431 The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian(). 2432 2433 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant() 2434 @*/ 2435 PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec X,Vec F,void *ctx) 2436 { 2437 PetscErrorCode ierr; 2438 Mat Arhs,Brhs; 2439 MatStructure flg2; 2440 2441 PetscFunctionBegin; 2442 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 2443 ierr = TSComputeRHSJacobian(ts,t,X,&Arhs,&Brhs,&flg2);CHKERRQ(ierr); 2444 ierr = MatMult(Arhs,X,F);CHKERRQ(ierr); 2445 PetscFunctionReturn(0); 2446 } 2447 2448 #undef __FUNCT__ 2449 #define __FUNCT__ "TSComputeRHSJacobianConstant" 2450 /*@C 2451 TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent. 2452 2453 Collective on TS 2454 2455 Input Arguments: 2456 + ts - time stepping context 2457 . t - time at which to evaluate 2458 . X - state at which to evaluate 2459 - ctx - context 2460 2461 Output Arguments: 2462 + A - pointer to operator 2463 . B - pointer to preconditioning matrix 2464 - flg - matrix structure flag 2465 2466 Level: intermediate 2467 2468 Notes: 2469 This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems. 2470 2471 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear() 2472 @*/ 2473 PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec X,Mat *A,Mat *B,MatStructure *flg,void *ctx) 2474 { 2475 2476 PetscFunctionBegin; 2477 *flg = SAME_PRECONDITIONER; 2478 PetscFunctionReturn(0); 2479 } 2480 2481 #undef __FUNCT__ 2482 #define __FUNCT__ "TSComputeIFunctionLinear" 2483 /*@C 2484 TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only 2485 2486 Collective on TS 2487 2488 Input Arguments: 2489 + ts - time stepping context 2490 . t - time at which to evaluate 2491 . X - state at which to evaluate 2492 . Xdot - time derivative of state vector 2493 - ctx - context 2494 2495 Output Arguments: 2496 . F - left hand side 2497 2498 Level: intermediate 2499 2500 Notes: 2501 The assumption here is that the left hand side is of the form A*Xdot (and not A*Xdot + B*X). For other cases, the 2502 user is required to write their own TSComputeIFunction. 2503 This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems. 2504 The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian(). 2505 2506 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant() 2507 @*/ 2508 PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx) 2509 { 2510 PetscErrorCode ierr; 2511 Mat A,B; 2512 MatStructure flg2; 2513 2514 PetscFunctionBegin; 2515 ierr = TSGetIJacobian(ts,&A,&B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2516 ierr = TSComputeIJacobian(ts,t,X,Xdot,1.0,&A,&B,&flg2,PETSC_TRUE);CHKERRQ(ierr); 2517 ierr = MatMult(A,Xdot,F);CHKERRQ(ierr); 2518 PetscFunctionReturn(0); 2519 } 2520 2521 #undef __FUNCT__ 2522 #define __FUNCT__ "TSComputeIJacobianConstant" 2523 /*@C 2524 TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent. 2525 2526 Collective on TS 2527 2528 Input Arguments: 2529 + ts - time stepping context 2530 . t - time at which to evaluate 2531 . X - state at which to evaluate 2532 . Xdot - time derivative of state vector 2533 . shift - shift to apply 2534 - ctx - context 2535 2536 Output Arguments: 2537 + A - pointer to operator 2538 . B - pointer to preconditioning matrix 2539 - flg - matrix structure flag 2540 2541 Level: intermediate 2542 2543 Notes: 2544 This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems. 2545 2546 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear() 2547 @*/ 2548 PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,void *ctx) 2549 { 2550 2551 PetscFunctionBegin; 2552 *flg = SAME_PRECONDITIONER; 2553 PetscFunctionReturn(0); 2554 } 2555 2556 2557 #undef __FUNCT__ 2558 #define __FUNCT__ "TSGetConvergedReason" 2559 /*@ 2560 TSGetConvergedReason - Gets the reason the TS iteration was stopped. 2561 2562 Not Collective 2563 2564 Input Parameter: 2565 . ts - the TS context 2566 2567 Output Parameter: 2568 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 2569 manual pages for the individual convergence tests for complete lists 2570 2571 Level: intermediate 2572 2573 Notes: 2574 Can only be called after the call to TSSolve() is complete. 2575 2576 .keywords: TS, nonlinear, set, convergence, test 2577 2578 .seealso: TSSetConvergenceTest(), TSConvergedReason 2579 @*/ 2580 PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason) 2581 { 2582 PetscFunctionBegin; 2583 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2584 PetscValidPointer(reason,2); 2585 *reason = ts->reason; 2586 PetscFunctionReturn(0); 2587 } 2588 2589 #if defined(PETSC_HAVE_MATLAB_ENGINE) 2590 #include <mex.h> 2591 2592 typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext; 2593 2594 #undef __FUNCT__ 2595 #define __FUNCT__ "TSComputeFunction_Matlab" 2596 /* 2597 TSComputeFunction_Matlab - Calls the function that has been set with 2598 TSSetFunctionMatlab(). 2599 2600 Collective on TS 2601 2602 Input Parameters: 2603 + snes - the TS context 2604 - x - input vector 2605 2606 Output Parameter: 2607 . y - function vector, as set by TSSetFunction() 2608 2609 Notes: 2610 TSComputeFunction() is typically used within nonlinear solvers 2611 implementations, so most users would not generally call this routine 2612 themselves. 2613 2614 Level: developer 2615 2616 .keywords: TS, nonlinear, compute, function 2617 2618 .seealso: TSSetFunction(), TSGetFunction() 2619 */ 2620 PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec x,Vec xdot,Vec y, void *ctx) 2621 { 2622 PetscErrorCode ierr; 2623 TSMatlabContext *sctx = (TSMatlabContext *)ctx; 2624 int nlhs = 1,nrhs = 7; 2625 mxArray *plhs[1],*prhs[7]; 2626 long long int lx = 0,lxdot = 0,ly = 0,ls = 0; 2627 2628 PetscFunctionBegin; 2629 PetscValidHeaderSpecific(snes,TS_CLASSID,1); 2630 PetscValidHeaderSpecific(x,VEC_CLASSID,3); 2631 PetscValidHeaderSpecific(xdot,VEC_CLASSID,4); 2632 PetscValidHeaderSpecific(y,VEC_CLASSID,5); 2633 PetscCheckSameComm(snes,1,x,3); 2634 PetscCheckSameComm(snes,1,y,5); 2635 2636 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 2637 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 2638 ierr = PetscMemcpy(&lxdot,&xdot,sizeof(xdot));CHKERRQ(ierr); 2639 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 2640 prhs[0] = mxCreateDoubleScalar((double)ls); 2641 prhs[1] = mxCreateDoubleScalar(time); 2642 prhs[2] = mxCreateDoubleScalar((double)lx); 2643 prhs[3] = mxCreateDoubleScalar((double)lxdot); 2644 prhs[4] = mxCreateDoubleScalar((double)ly); 2645 prhs[5] = mxCreateString(sctx->funcname); 2646 prhs[6] = sctx->ctx; 2647 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");CHKERRQ(ierr); 2648 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 2649 mxDestroyArray(prhs[0]); 2650 mxDestroyArray(prhs[1]); 2651 mxDestroyArray(prhs[2]); 2652 mxDestroyArray(prhs[3]); 2653 mxDestroyArray(prhs[4]); 2654 mxDestroyArray(prhs[5]); 2655 mxDestroyArray(plhs[0]); 2656 PetscFunctionReturn(0); 2657 } 2658 2659 2660 #undef __FUNCT__ 2661 #define __FUNCT__ "TSSetFunctionMatlab" 2662 /* 2663 TSSetFunctionMatlab - Sets the function evaluation routine and function 2664 vector for use by the TS routines in solving ODEs 2665 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 2666 2667 Logically Collective on TS 2668 2669 Input Parameters: 2670 + ts - the TS context 2671 - func - function evaluation routine 2672 2673 Calling sequence of func: 2674 $ func (TS ts,PetscReal time,Vec x,Vec xdot,Vec f,void *ctx); 2675 2676 Level: beginner 2677 2678 .keywords: TS, nonlinear, set, function 2679 2680 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 2681 */ 2682 PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx) 2683 { 2684 PetscErrorCode ierr; 2685 TSMatlabContext *sctx; 2686 2687 PetscFunctionBegin; 2688 /* currently sctx is memory bleed */ 2689 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 2690 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 2691 /* 2692 This should work, but it doesn't 2693 sctx->ctx = ctx; 2694 mexMakeArrayPersistent(sctx->ctx); 2695 */ 2696 sctx->ctx = mxDuplicateArray(ctx); 2697 ierr = TSSetIFunction(ts,PETSC_NULL,TSComputeFunction_Matlab,sctx);CHKERRQ(ierr); 2698 PetscFunctionReturn(0); 2699 } 2700 2701 #undef __FUNCT__ 2702 #define __FUNCT__ "TSComputeJacobian_Matlab" 2703 /* 2704 TSComputeJacobian_Matlab - Calls the function that has been set with 2705 TSSetJacobianMatlab(). 2706 2707 Collective on TS 2708 2709 Input Parameters: 2710 + ts - the TS context 2711 . x - input vector 2712 . A, B - the matrices 2713 - ctx - user context 2714 2715 Output Parameter: 2716 . flag - structure of the matrix 2717 2718 Level: developer 2719 2720 .keywords: TS, nonlinear, compute, function 2721 2722 .seealso: TSSetFunction(), TSGetFunction() 2723 @*/ 2724 PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec x,Vec xdot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag, void *ctx) 2725 { 2726 PetscErrorCode ierr; 2727 TSMatlabContext *sctx = (TSMatlabContext *)ctx; 2728 int nlhs = 2,nrhs = 9; 2729 mxArray *plhs[2],*prhs[9]; 2730 long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0; 2731 2732 PetscFunctionBegin; 2733 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2734 PetscValidHeaderSpecific(x,VEC_CLASSID,3); 2735 2736 /* call Matlab function in ctx with arguments x and y */ 2737 2738 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 2739 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 2740 ierr = PetscMemcpy(&lxdot,&xdot,sizeof(x));CHKERRQ(ierr); 2741 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 2742 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 2743 prhs[0] = mxCreateDoubleScalar((double)ls); 2744 prhs[1] = mxCreateDoubleScalar((double)time); 2745 prhs[2] = mxCreateDoubleScalar((double)lx); 2746 prhs[3] = mxCreateDoubleScalar((double)lxdot); 2747 prhs[4] = mxCreateDoubleScalar((double)shift); 2748 prhs[5] = mxCreateDoubleScalar((double)lA); 2749 prhs[6] = mxCreateDoubleScalar((double)lB); 2750 prhs[7] = mxCreateString(sctx->funcname); 2751 prhs[8] = sctx->ctx; 2752 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");CHKERRQ(ierr); 2753 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 2754 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 2755 mxDestroyArray(prhs[0]); 2756 mxDestroyArray(prhs[1]); 2757 mxDestroyArray(prhs[2]); 2758 mxDestroyArray(prhs[3]); 2759 mxDestroyArray(prhs[4]); 2760 mxDestroyArray(prhs[5]); 2761 mxDestroyArray(prhs[6]); 2762 mxDestroyArray(prhs[7]); 2763 mxDestroyArray(plhs[0]); 2764 mxDestroyArray(plhs[1]); 2765 PetscFunctionReturn(0); 2766 } 2767 2768 2769 #undef __FUNCT__ 2770 #define __FUNCT__ "TSSetJacobianMatlab" 2771 /* 2772 TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 2773 vector for use by the TS routines in solving ODEs from MATLAB. Here the function is a string containing the name of a MATLAB function 2774 2775 Logically Collective on TS 2776 2777 Input Parameters: 2778 + ts - the TS context 2779 . A,B - Jacobian matrices 2780 . func - function evaluation routine 2781 - ctx - user context 2782 2783 Calling sequence of func: 2784 $ flag = func (TS ts,PetscReal time,Vec x,Vec xdot,Mat A,Mat B,void *ctx); 2785 2786 2787 Level: developer 2788 2789 .keywords: TS, nonlinear, set, function 2790 2791 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 2792 */ 2793 PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx) 2794 { 2795 PetscErrorCode ierr; 2796 TSMatlabContext *sctx; 2797 2798 PetscFunctionBegin; 2799 /* currently sctx is memory bleed */ 2800 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 2801 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 2802 /* 2803 This should work, but it doesn't 2804 sctx->ctx = ctx; 2805 mexMakeArrayPersistent(sctx->ctx); 2806 */ 2807 sctx->ctx = mxDuplicateArray(ctx); 2808 ierr = TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 2809 PetscFunctionReturn(0); 2810 } 2811 2812 #undef __FUNCT__ 2813 #define __FUNCT__ "TSMonitor_Matlab" 2814 /* 2815 TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab(). 2816 2817 Collective on TS 2818 2819 .seealso: TSSetFunction(), TSGetFunction() 2820 @*/ 2821 PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec x, void *ctx) 2822 { 2823 PetscErrorCode ierr; 2824 TSMatlabContext *sctx = (TSMatlabContext *)ctx; 2825 int nlhs = 1,nrhs = 6; 2826 mxArray *plhs[1],*prhs[6]; 2827 long long int lx = 0,ls = 0; 2828 2829 PetscFunctionBegin; 2830 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2831 PetscValidHeaderSpecific(x,VEC_CLASSID,4); 2832 2833 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 2834 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 2835 prhs[0] = mxCreateDoubleScalar((double)ls); 2836 prhs[1] = mxCreateDoubleScalar((double)it); 2837 prhs[2] = mxCreateDoubleScalar((double)time); 2838 prhs[3] = mxCreateDoubleScalar((double)lx); 2839 prhs[4] = mxCreateString(sctx->funcname); 2840 prhs[5] = sctx->ctx; 2841 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");CHKERRQ(ierr); 2842 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 2843 mxDestroyArray(prhs[0]); 2844 mxDestroyArray(prhs[1]); 2845 mxDestroyArray(prhs[2]); 2846 mxDestroyArray(prhs[3]); 2847 mxDestroyArray(prhs[4]); 2848 mxDestroyArray(plhs[0]); 2849 PetscFunctionReturn(0); 2850 } 2851 2852 2853 #undef __FUNCT__ 2854 #define __FUNCT__ "TSMonitorSetMatlab" 2855 /* 2856 TSMonitorSetMatlab - Sets the monitor function from Matlab 2857 2858 Level: developer 2859 2860 .keywords: TS, nonlinear, set, function 2861 2862 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 2863 */ 2864 PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx) 2865 { 2866 PetscErrorCode ierr; 2867 TSMatlabContext *sctx; 2868 2869 PetscFunctionBegin; 2870 /* currently sctx is memory bleed */ 2871 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 2872 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 2873 /* 2874 This should work, but it doesn't 2875 sctx->ctx = ctx; 2876 mexMakeArrayPersistent(sctx->ctx); 2877 */ 2878 sctx->ctx = mxDuplicateArray(ctx); 2879 ierr = TSMonitorSet(ts,TSMonitor_Matlab,sctx,PETSC_NULL);CHKERRQ(ierr); 2880 PetscFunctionReturn(0); 2881 } 2882 2883 #endif 2884