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