1 2 #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/ 3 #include <petscdmshell.h> 4 5 /* Logging support */ 6 PetscClassId TS_CLASSID; 7 PetscLogEvent TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval; 8 9 #undef __FUNCT__ 10 #define __FUNCT__ "TSSetTypeFromOptions" 11 /* 12 TSSetTypeFromOptions - Sets the type of ts from user options. 13 14 Collective on TS 15 16 Input Parameter: 17 . ts - The ts 18 19 Level: intermediate 20 21 .keywords: TS, set, options, database, type 22 .seealso: TSSetFromOptions(), TSSetType() 23 */ 24 static PetscErrorCode TSSetTypeFromOptions(TS ts) 25 { 26 PetscBool opt; 27 const char *defaultType; 28 char typeName[256]; 29 PetscErrorCode ierr; 30 31 PetscFunctionBegin; 32 if (((PetscObject)ts)->type_name) { 33 defaultType = ((PetscObject)ts)->type_name; 34 } else { 35 defaultType = TSEULER; 36 } 37 38 if (!TSRegisterAllCalled) {ierr = TSRegisterAll(PETSC_NULL);CHKERRQ(ierr);} 39 ierr = PetscOptionsList("-ts_type", "TS method"," TSSetType", TSList, defaultType, typeName, 256, &opt);CHKERRQ(ierr); 40 if (opt) { 41 ierr = TSSetType(ts, typeName);CHKERRQ(ierr); 42 } else { 43 ierr = TSSetType(ts, defaultType);CHKERRQ(ierr); 44 } 45 PetscFunctionReturn(0); 46 } 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "TSSetFromOptions" 50 /*@ 51 TSSetFromOptions - Sets various TS parameters from user options. 52 53 Collective on TS 54 55 Input Parameter: 56 . ts - the TS context obtained from TSCreate() 57 58 Options Database Keys: 59 + -ts_type <type> - TSEULER, TSBEULER, TSSUNDIALS, TSPSEUDO, TSCN, TSRK, TSTHETA, TSGL, TSSSP 60 . -ts_max_steps maxsteps - maximum number of time-steps to take 61 . -ts_final_time time - maximum time to compute to 62 . -ts_dt dt - initial time step 63 . -ts_monitor - print information at each timestep 64 - -ts_monitor_draw - plot information at each timestep 65 66 Level: beginner 67 68 .keywords: TS, timestep, set, options, database 69 70 .seealso: TSGetType() 71 @*/ 72 PetscErrorCode TSSetFromOptions(TS ts) 73 { 74 PetscBool opt,flg; 75 PetscErrorCode ierr; 76 PetscViewer monviewer; 77 char monfilename[PETSC_MAX_PATH_LEN]; 78 SNES snes; 79 TSAdapt adapt; 80 PetscReal time_step; 81 82 PetscFunctionBegin; 83 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 84 ierr = PetscObjectOptionsBegin((PetscObject)ts);CHKERRQ(ierr); 85 /* Handle TS type options */ 86 ierr = TSSetTypeFromOptions(ts);CHKERRQ(ierr); 87 88 /* Handle generic TS options */ 89 ierr = PetscOptionsInt("-ts_max_steps","Maximum number of time steps","TSSetDuration",ts->max_steps,&ts->max_steps,PETSC_NULL);CHKERRQ(ierr); 90 ierr = PetscOptionsReal("-ts_final_time","Time to run to","TSSetDuration",ts->max_time,&ts->max_time,PETSC_NULL);CHKERRQ(ierr); 91 ierr = PetscOptionsReal("-ts_init_time","Initial time","TSSetTime",ts->ptime,&ts->ptime,PETSC_NULL);CHKERRQ(ierr); 92 ierr = PetscOptionsReal("-ts_dt","Initial time step","TSSetTimeStep",ts->time_step,&time_step,&flg);CHKERRQ(ierr); 93 if (flg) { 94 ierr = TSSetTimeStep(ts,time_step);CHKERRQ(ierr); 95 } 96 opt = ts->exact_final_time == PETSC_DECIDE ? PETSC_FALSE : (PetscBool)ts->exact_final_time; 97 ierr = PetscOptionsBool("-ts_exact_final_time","Interpolate output to stop exactly at the final time","TSSetExactFinalTime",opt,&opt,&flg);CHKERRQ(ierr); 98 if (flg) {ierr = TSSetExactFinalTime(ts,opt);CHKERRQ(ierr);} 99 ierr = PetscOptionsInt("-ts_max_snes_failures","Maximum number of nonlinear solve failures","TSSetMaxSNESFailures",ts->max_snes_failures,&ts->max_snes_failures,PETSC_NULL);CHKERRQ(ierr); 100 ierr = PetscOptionsInt("-ts_max_reject","Maximum number of step rejections before step fails","TSSetMaxStepRejections",ts->max_reject,&ts->max_reject,PETSC_NULL);CHKERRQ(ierr); 101 ierr = PetscOptionsBool("-ts_error_if_step_fails","Error if no step succeeds","TSSetErrorIfStepFails",ts->errorifstepfailed,&ts->errorifstepfailed,PETSC_NULL);CHKERRQ(ierr); 102 ierr = PetscOptionsReal("-ts_rtol","Relative tolerance for local truncation error","TSSetTolerances",ts->rtol,&ts->rtol,PETSC_NULL);CHKERRQ(ierr); 103 ierr = PetscOptionsReal("-ts_atol","Absolute tolerance for local truncation error","TSSetTolerances",ts->atol,&ts->atol,PETSC_NULL);CHKERRQ(ierr); 104 105 /* Monitor options */ 106 ierr = PetscOptionsString("-ts_monitor","Monitor timestep size","TSMonitorDefault","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 107 if (flg) { 108 ierr = PetscViewerASCIIOpen(((PetscObject)ts)->comm,monfilename,&monviewer);CHKERRQ(ierr); 109 ierr = TSMonitorSet(ts,TSMonitorDefault,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);CHKERRQ(ierr); 110 } 111 ierr = PetscOptionsString("-ts_monitor_python","Use Python function","TSMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 112 if (flg) {ierr = PetscPythonMonitorSet((PetscObject)ts,monfilename);CHKERRQ(ierr);} 113 114 opt = PETSC_FALSE; 115 ierr = PetscOptionsBool("-ts_monitor_draw","Monitor timestep size graphically","TSMonitorLG",opt,&opt,PETSC_NULL);CHKERRQ(ierr); 116 if (opt) { 117 PetscDrawLG lg; 118 119 ierr = TSMonitorLGCreate(((PetscObject)ts)->comm,0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,&lg); 120 ierr = TSMonitorSet(ts,TSMonitorLG,lg,(PetscErrorCode (*)(void**))TSMonitorLGDestroy);CHKERRQ(ierr); 121 } 122 opt = PETSC_FALSE; 123 ierr = PetscOptionsBool("-ts_monitor_solutionode","Monitor solution graphically","TSMonitorSolutionODE",opt,&opt,PETSC_NULL);CHKERRQ(ierr); 124 if (opt) { 125 PetscDrawLG lg; 126 127 ierr = TSMonitorSolutionODECreate(((PetscObject)ts)->comm,3,0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,&lg); 128 ierr = TSMonitorSet(ts,TSMonitorSolutionODE,lg,(PetscErrorCode (*)(void**))TSMonitorSolutionODEDestroy);CHKERRQ(ierr); 129 } 130 opt = PETSC_FALSE; 131 ierr = PetscOptionsBool("-ts_monitor_errorode","Monitor error graphically","TSMonitorErrorODE",opt,&opt,PETSC_NULL);CHKERRQ(ierr); 132 if (opt) { 133 PetscDrawLG lg; 134 135 ierr = TSMonitorErrorODECreate(((PetscObject)ts)->comm,3,0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,&lg); 136 ierr = TSMonitorSet(ts,TSMonitorErrorODE,lg,(PetscErrorCode (*)(void**))TSMonitorErrorODEDestroy);CHKERRQ(ierr); 137 } 138 opt = PETSC_FALSE; 139 ierr = PetscOptionsBool("-ts_monitor_solution","Monitor solution graphically","TSMonitorSolution",opt,&opt,PETSC_NULL);CHKERRQ(ierr); 140 if (opt) { 141 void *ctx; 142 PetscViewer viewer; 143 144 ierr = PetscViewerDrawOpen(((PetscObject)ts)->comm,0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,&viewer); 145 ierr = TSMonitorSolutionCreate(ts,viewer,&ctx);CHKERRQ(ierr); 146 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 147 ierr = TSMonitorSet(ts,TSMonitorSolution,ctx,TSMonitorSolutionDestroy);CHKERRQ(ierr); 148 } 149 opt = PETSC_FALSE; 150 ierr = PetscOptionsString("-ts_monitor_solution_binary","Save each solution to a binary file","TSMonitorSolutionBinary",0,monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 151 if (flg) { 152 PetscViewer ctx; 153 if (monfilename[0]) { 154 ierr = PetscViewerBinaryOpen(((PetscObject)ts)->comm,monfilename,FILE_MODE_WRITE,&ctx);CHKERRQ(ierr); 155 } else { 156 ctx = PETSC_VIEWER_BINARY_(((PetscObject)ts)->comm); 157 } 158 ierr = TSMonitorSet(ts,TSMonitorSolutionBinary,ctx,(PetscErrorCode (*)(void**))PetscViewerDestroy);CHKERRQ(ierr); 159 } 160 opt = PETSC_FALSE; 161 ierr = PetscOptionsString("-ts_monitor_solution_vtk","Save each time step to a binary file, use filename-%%03D.vts","TSMonitorSolutionVTK",0,monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 162 if (flg) { 163 const char *ptr,*ptr2; 164 char *filetemplate; 165 if (!monfilename[0]) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts"); 166 /* Do some cursory validation of the input. */ 167 ierr = PetscStrstr(monfilename,"%",(char**)&ptr);CHKERRQ(ierr); 168 if (!ptr) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts"); 169 for (ptr++ ; ptr && *ptr; ptr++) { 170 ierr = PetscStrchr("DdiouxX",*ptr,(char**)&ptr2);CHKERRQ(ierr); 171 if (!ptr2 && (*ptr < '0' || '9' < *ptr)) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"Invalid file template argument to -ts_monitor_solution_vtk, should look like filename-%%03D.vts"); 172 if (ptr2) break; 173 } 174 ierr = PetscStrallocpy(monfilename,&filetemplate);CHKERRQ(ierr); 175 ierr = TSMonitorSet(ts,TSMonitorSolutionVTK,filetemplate,(PetscErrorCode (*)(void**))TSMonitorSolutionVTKDestroy);CHKERRQ(ierr); 176 } 177 178 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 179 ierr = TSAdaptSetFromOptions(adapt);CHKERRQ(ierr); 180 181 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 182 if (ts->problem_type == TS_LINEAR) {ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr);} 183 184 /* Handle specific TS options */ 185 if (ts->ops->setfromoptions) { 186 ierr = (*ts->ops->setfromoptions)(ts);CHKERRQ(ierr); 187 } 188 189 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 190 ierr = PetscObjectProcessOptionsHandlers((PetscObject)ts);CHKERRQ(ierr); 191 ierr = PetscOptionsEnd();CHKERRQ(ierr); 192 PetscFunctionReturn(0); 193 } 194 195 #undef __FUNCT__ 196 #undef __FUNCT__ 197 #define __FUNCT__ "TSComputeRHSJacobian" 198 /*@ 199 TSComputeRHSJacobian - Computes the Jacobian matrix that has been 200 set with TSSetRHSJacobian(). 201 202 Collective on TS and Vec 203 204 Input Parameters: 205 + ts - the TS context 206 . t - current timestep 207 - x - input vector 208 209 Output Parameters: 210 + A - Jacobian matrix 211 . B - optional preconditioning matrix 212 - flag - flag indicating matrix structure 213 214 Notes: 215 Most users should not need to explicitly call this routine, as it 216 is used internally within the nonlinear solvers. 217 218 See KSPSetOperators() for important information about setting the 219 flag parameter. 220 221 Level: developer 222 223 .keywords: SNES, compute, Jacobian, matrix 224 225 .seealso: TSSetRHSJacobian(), KSPSetOperators() 226 @*/ 227 PetscErrorCode TSComputeRHSJacobian(TS ts,PetscReal t,Vec X,Mat *A,Mat *B,MatStructure *flg) 228 { 229 PetscErrorCode ierr; 230 PetscInt Xstate; 231 DM dm; 232 TSDM tsdm; 233 TSRHSJacobian rhsjacobianfunc; 234 void *ctx; 235 TSIJacobian ijacobianfunc; 236 237 PetscFunctionBegin; 238 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 239 PetscValidHeaderSpecific(X,VEC_CLASSID,3); 240 PetscCheckSameComm(ts,1,X,3); 241 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 242 ierr = DMTSGetContext(dm,&tsdm);CHKERRQ(ierr); 243 ierr = DMTSGetRHSJacobian(dm,&rhsjacobianfunc,&ctx);CHKERRQ(ierr); 244 ierr = DMTSGetIJacobian(dm,&ijacobianfunc,PETSC_NULL);CHKERRQ(ierr); 245 ierr = PetscObjectStateQuery((PetscObject)X,&Xstate);CHKERRQ(ierr); 246 if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.X == X && ts->rhsjacobian.Xstate == Xstate))) { 247 *flg = ts->rhsjacobian.mstructure; 248 PetscFunctionReturn(0); 249 } 250 251 if (!rhsjacobianfunc && !ijacobianfunc) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()"); 252 253 if (rhsjacobianfunc) { 254 ierr = PetscLogEventBegin(TS_JacobianEval,ts,X,*A,*B);CHKERRQ(ierr); 255 *flg = DIFFERENT_NONZERO_PATTERN; 256 PetscStackPush("TS user Jacobian function"); 257 ierr = (*rhsjacobianfunc)(ts,t,X,A,B,flg,ctx);CHKERRQ(ierr); 258 PetscStackPop; 259 ierr = PetscLogEventEnd(TS_JacobianEval,ts,X,*A,*B);CHKERRQ(ierr); 260 /* make sure user returned a correct Jacobian and preconditioner */ 261 PetscValidHeaderSpecific(*A,MAT_CLASSID,4); 262 PetscValidHeaderSpecific(*B,MAT_CLASSID,5); 263 } else { 264 ierr = MatZeroEntries(*A);CHKERRQ(ierr); 265 if (*A != *B) {ierr = MatZeroEntries(*B);CHKERRQ(ierr);} 266 *flg = SAME_NONZERO_PATTERN; 267 } 268 ts->rhsjacobian.time = t; 269 ts->rhsjacobian.X = X; 270 ierr = PetscObjectStateQuery((PetscObject)X,&ts->rhsjacobian.Xstate);CHKERRQ(ierr); 271 ts->rhsjacobian.mstructure = *flg; 272 PetscFunctionReturn(0); 273 } 274 275 #undef __FUNCT__ 276 #define __FUNCT__ "TSComputeRHSFunction" 277 /*@ 278 TSComputeRHSFunction - Evaluates the right-hand-side function. 279 280 Collective on TS and Vec 281 282 Input Parameters: 283 + ts - the TS context 284 . t - current time 285 - x - state vector 286 287 Output Parameter: 288 . y - right hand side 289 290 Note: 291 Most users should not need to explicitly call this routine, as it 292 is used internally within the nonlinear solvers. 293 294 Level: developer 295 296 .keywords: TS, compute 297 298 .seealso: TSSetRHSFunction(), TSComputeIFunction() 299 @*/ 300 PetscErrorCode TSComputeRHSFunction(TS ts,PetscReal t,Vec x,Vec y) 301 { 302 PetscErrorCode ierr; 303 304 TSRHSFunction rhsfunction; 305 TSIFunction ifunction; 306 void *ctx; 307 DM dm; 308 309 PetscFunctionBegin; 310 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 311 PetscValidHeaderSpecific(x,VEC_CLASSID,3); 312 PetscValidHeaderSpecific(y,VEC_CLASSID,4); 313 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 314 ierr = DMTSGetRHSFunction(dm,&rhsfunction,&ctx);CHKERRQ(ierr); 315 ierr = DMTSGetIFunction(dm,&ifunction,PETSC_NULL);CHKERRQ(ierr); 316 317 if (!rhsfunction && !ifunction) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()"); 318 319 ierr = PetscLogEventBegin(TS_FunctionEval,ts,x,y,0);CHKERRQ(ierr); 320 if (rhsfunction) { 321 PetscStackPush("TS user right-hand-side function"); 322 ierr = (*rhsfunction)(ts,t,x,y,ctx);CHKERRQ(ierr); 323 PetscStackPop; 324 } else { 325 ierr = VecZeroEntries(y);CHKERRQ(ierr); 326 } 327 328 ierr = PetscLogEventEnd(TS_FunctionEval,ts,x,y,0);CHKERRQ(ierr); 329 PetscFunctionReturn(0); 330 } 331 332 #undef __FUNCT__ 333 #define __FUNCT__ "TSComputeSolutionFunction" 334 /*@ 335 TSComputeSolutionFunction - Evaluates the solution function. 336 337 Collective on TS and Vec 338 339 Input Parameters: 340 + ts - the TS context 341 - t - current time 342 343 Output Parameter: 344 . y - right hand side 345 346 Note: 347 Most users should not need to explicitly call this routine, as it 348 is used internally within the nonlinear solvers. 349 350 Level: developer 351 352 .keywords: TS, compute 353 354 .seealso: TSSetSolutionFunction(), TSSetRHSFunction(), TSComputeIFunction() 355 @*/ 356 PetscErrorCode TSComputeSolutionFunction(TS ts,PetscReal t,Vec x) 357 { 358 PetscErrorCode ierr; 359 TSSolutionFunction solutionfunction; 360 void *ctx; 361 DM dm; 362 363 PetscFunctionBegin; 364 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 365 PetscValidHeaderSpecific(x,VEC_CLASSID,3); 366 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 367 ierr = DMTSGetSolutionFunction(dm,&solutionfunction,&ctx);CHKERRQ(ierr); 368 369 if (solutionfunction) { 370 PetscStackPush("TS user right-hand-side function"); 371 ierr = (*solutionfunction)(ts,t,x,ctx);CHKERRQ(ierr); 372 PetscStackPop; 373 } 374 PetscFunctionReturn(0); 375 } 376 377 #undef __FUNCT__ 378 #define __FUNCT__ "TSGetRHSVec_Private" 379 static PetscErrorCode TSGetRHSVec_Private(TS ts,Vec *Frhs) 380 { 381 Vec F; 382 PetscErrorCode ierr; 383 384 PetscFunctionBegin; 385 *Frhs = PETSC_NULL; 386 ierr = TSGetIFunction(ts,&F,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 387 if (!ts->Frhs) { 388 ierr = VecDuplicate(F,&ts->Frhs);CHKERRQ(ierr); 389 } 390 *Frhs = ts->Frhs; 391 PetscFunctionReturn(0); 392 } 393 394 #undef __FUNCT__ 395 #define __FUNCT__ "TSGetRHSMats_Private" 396 static PetscErrorCode TSGetRHSMats_Private(TS ts,Mat *Arhs,Mat *Brhs) 397 { 398 Mat A,B; 399 PetscErrorCode ierr; 400 401 PetscFunctionBegin; 402 ierr = TSGetIJacobian(ts,&A,&B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 403 if (Arhs) { 404 if (!ts->Arhs) { 405 ierr = MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&ts->Arhs);CHKERRQ(ierr); 406 } 407 *Arhs = ts->Arhs; 408 } 409 if (Brhs) { 410 if (!ts->Brhs) { 411 ierr = MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&ts->Brhs);CHKERRQ(ierr); 412 } 413 *Brhs = ts->Brhs; 414 } 415 PetscFunctionReturn(0); 416 } 417 418 #undef __FUNCT__ 419 #define __FUNCT__ "TSComputeIFunction" 420 /*@ 421 TSComputeIFunction - Evaluates the DAE residual written in implicit form F(t,X,Xdot)=0 422 423 Collective on TS and Vec 424 425 Input Parameters: 426 + ts - the TS context 427 . t - current time 428 . X - state vector 429 . Xdot - time derivative of state vector 430 - imex - flag indicates if the method is IMEX so that the RHSFunction should be kept separate 431 432 Output Parameter: 433 . Y - right hand side 434 435 Note: 436 Most users should not need to explicitly call this routine, as it 437 is used internally within the nonlinear solvers. 438 439 If the user did did not write their equations in implicit form, this 440 function recasts them in implicit form. 441 442 Level: developer 443 444 .keywords: TS, compute 445 446 .seealso: TSSetIFunction(), TSComputeRHSFunction() 447 @*/ 448 PetscErrorCode TSComputeIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec Y,PetscBool imex) 449 { 450 PetscErrorCode ierr; 451 TSIFunction ifunction; 452 TSRHSFunction rhsfunction; 453 void *ctx; 454 DM dm; 455 456 PetscFunctionBegin; 457 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 458 PetscValidHeaderSpecific(X,VEC_CLASSID,3); 459 PetscValidHeaderSpecific(Xdot,VEC_CLASSID,4); 460 PetscValidHeaderSpecific(Y,VEC_CLASSID,5); 461 462 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 463 ierr = DMTSGetIFunction(dm,&ifunction,&ctx);CHKERRQ(ierr); 464 ierr = DMTSGetRHSFunction(dm,&rhsfunction,PETSC_NULL);CHKERRQ(ierr); 465 466 if (!rhsfunction && !ifunction) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()"); 467 468 ierr = PetscLogEventBegin(TS_FunctionEval,ts,X,Xdot,Y);CHKERRQ(ierr); 469 if (ifunction) { 470 PetscStackPush("TS user implicit function"); 471 ierr = (*ifunction)(ts,t,X,Xdot,Y,ctx);CHKERRQ(ierr); 472 PetscStackPop; 473 } 474 if (imex) { 475 if (!ifunction) { 476 ierr = VecCopy(Xdot,Y);CHKERRQ(ierr); 477 } 478 } else if (rhsfunction) { 479 if (ifunction) { 480 Vec Frhs; 481 ierr = TSGetRHSVec_Private(ts,&Frhs);CHKERRQ(ierr); 482 ierr = TSComputeRHSFunction(ts,t,X,Frhs);CHKERRQ(ierr); 483 ierr = VecAXPY(Y,-1,Frhs);CHKERRQ(ierr); 484 } else { 485 ierr = TSComputeRHSFunction(ts,t,X,Y);CHKERRQ(ierr); 486 ierr = VecAYPX(Y,-1,Xdot);CHKERRQ(ierr); 487 } 488 } 489 ierr = PetscLogEventEnd(TS_FunctionEval,ts,X,Xdot,Y);CHKERRQ(ierr); 490 PetscFunctionReturn(0); 491 } 492 493 #undef __FUNCT__ 494 #define __FUNCT__ "TSComputeIJacobian" 495 /*@ 496 TSComputeIJacobian - Evaluates the Jacobian of the DAE 497 498 Collective on TS and Vec 499 500 Input 501 Input Parameters: 502 + ts - the TS context 503 . t - current timestep 504 . X - state vector 505 . Xdot - time derivative of state vector 506 . shift - shift to apply, see note below 507 - imex - flag indicates if the method is IMEX so that the RHSJacobian should be kept separate 508 509 Output Parameters: 510 + A - Jacobian matrix 511 . B - optional preconditioning matrix 512 - flag - flag indicating matrix structure 513 514 Notes: 515 If F(t,X,Xdot)=0 is the DAE, the required Jacobian is 516 517 dF/dX + shift*dF/dXdot 518 519 Most users should not need to explicitly call this routine, as it 520 is used internally within the nonlinear solvers. 521 522 Level: developer 523 524 .keywords: TS, compute, Jacobian, matrix 525 526 .seealso: TSSetIJacobian() 527 @*/ 528 PetscErrorCode TSComputeIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,PetscBool imex) 529 { 530 PetscInt Xstate, Xdotstate; 531 PetscErrorCode ierr; 532 TSIJacobian ijacobian; 533 TSRHSJacobian rhsjacobian; 534 DM dm; 535 void *ctx; 536 537 PetscFunctionBegin; 538 539 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 540 PetscValidHeaderSpecific(X,VEC_CLASSID,3); 541 PetscValidHeaderSpecific(Xdot,VEC_CLASSID,4); 542 PetscValidPointer(A,6); 543 PetscValidHeaderSpecific(*A,MAT_CLASSID,6); 544 PetscValidPointer(B,7); 545 PetscValidHeaderSpecific(*B,MAT_CLASSID,7); 546 PetscValidPointer(flg,8); 547 548 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 549 ierr = DMTSGetIJacobian(dm,&ijacobian,&ctx);CHKERRQ(ierr); 550 ierr = DMTSGetRHSJacobian(dm,&rhsjacobian,PETSC_NULL);CHKERRQ(ierr); 551 552 ierr = PetscObjectStateQuery((PetscObject)X,&Xstate);CHKERRQ(ierr); 553 ierr = PetscObjectStateQuery((PetscObject)Xdot,&Xdotstate);CHKERRQ(ierr); 554 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))) { 555 *flg = ts->ijacobian.mstructure; 556 ierr = MatScale(*A, shift / ts->ijacobian.shift);CHKERRQ(ierr); 557 PetscFunctionReturn(0); 558 } 559 560 if (!rhsjacobian && !ijacobian) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()"); 561 562 *flg = SAME_NONZERO_PATTERN; /* In case we're solving a linear problem in which case it wouldn't get initialized below. */ 563 ierr = PetscLogEventBegin(TS_JacobianEval,ts,X,*A,*B);CHKERRQ(ierr); 564 if (ijacobian) { 565 *flg = DIFFERENT_NONZERO_PATTERN; 566 PetscStackPush("TS user implicit Jacobian"); 567 ierr = (*ijacobian)(ts,t,X,Xdot,shift,A,B,flg,ctx);CHKERRQ(ierr); 568 PetscStackPop; 569 /* make sure user returned a correct Jacobian and preconditioner */ 570 PetscValidHeaderSpecific(*A,MAT_CLASSID,4); 571 PetscValidHeaderSpecific(*B,MAT_CLASSID,5); 572 } 573 if (imex) { 574 if (!ijacobian) { /* system was written as Xdot = F(t,X) */ 575 ierr = MatZeroEntries(*A);CHKERRQ(ierr); 576 ierr = MatShift(*A,shift);CHKERRQ(ierr); 577 if (*A != *B) { 578 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 579 ierr = MatShift(*B,shift);CHKERRQ(ierr); 580 } 581 *flg = SAME_PRECONDITIONER; 582 } 583 } else { 584 if (!ijacobian) { 585 ierr = TSComputeRHSJacobian(ts,t,X,A,B,flg);CHKERRQ(ierr); 586 ierr = MatScale(*A,-1);CHKERRQ(ierr); 587 ierr = MatShift(*A,shift);CHKERRQ(ierr); 588 if (*A != *B) { 589 ierr = MatScale(*B,-1);CHKERRQ(ierr); 590 ierr = MatShift(*B,shift);CHKERRQ(ierr); 591 } 592 } else if (rhsjacobian) { 593 Mat Arhs,Brhs; 594 MatStructure axpy,flg2 = DIFFERENT_NONZERO_PATTERN; 595 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 596 ierr = TSComputeRHSJacobian(ts,t,X,&Arhs,&Brhs,&flg2);CHKERRQ(ierr); 597 axpy = (*flg == flg2) ? SAME_NONZERO_PATTERN : DIFFERENT_NONZERO_PATTERN; 598 ierr = MatAXPY(*A,-1,Arhs,axpy);CHKERRQ(ierr); 599 if (*A != *B) { 600 ierr = MatAXPY(*B,-1,Brhs,axpy);CHKERRQ(ierr); 601 } 602 *flg = PetscMin(*flg,flg2); 603 } 604 } 605 606 ts->ijacobian.time = t; 607 ts->ijacobian.X = X; 608 ts->ijacobian.Xdot = Xdot; 609 ierr = PetscObjectStateQuery((PetscObject)X,&ts->ijacobian.Xstate);CHKERRQ(ierr); 610 ierr = PetscObjectStateQuery((PetscObject)Xdot,&ts->ijacobian.Xdotstate);CHKERRQ(ierr); 611 ts->ijacobian.shift = shift; 612 ts->ijacobian.imex = imex; 613 ts->ijacobian.mstructure = *flg; 614 ierr = PetscLogEventEnd(TS_JacobianEval,ts,X,*A,*B);CHKERRQ(ierr); 615 PetscFunctionReturn(0); 616 } 617 618 #undef __FUNCT__ 619 #define __FUNCT__ "TSSetRHSFunction" 620 /*@C 621 TSSetRHSFunction - Sets the routine for evaluating the function, 622 F(t,u), where U_t = F(t,u). 623 624 Logically Collective on TS 625 626 Input Parameters: 627 + ts - the TS context obtained from TSCreate() 628 . r - vector to put the computed right hand side (or PETSC_NULL to have it created) 629 . f - routine for evaluating the right-hand-side function 630 - ctx - [optional] user-defined context for private data for the 631 function evaluation routine (may be PETSC_NULL) 632 633 Calling sequence of func: 634 $ func (TS ts,PetscReal t,Vec u,Vec F,void *ctx); 635 636 + t - current timestep 637 . u - input vector 638 . F - function vector 639 - ctx - [optional] user-defined function context 640 641 Level: beginner 642 643 .keywords: TS, timestep, set, right-hand-side, function 644 645 .seealso: TSSetRHSJacobian(), TSSetIJacobian() 646 @*/ 647 PetscErrorCode TSSetRHSFunction(TS ts,Vec r,PetscErrorCode (*f)(TS,PetscReal,Vec,Vec,void*),void *ctx) 648 { 649 PetscErrorCode ierr; 650 SNES snes; 651 Vec ralloc = PETSC_NULL; 652 DM dm; 653 654 PetscFunctionBegin; 655 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 656 if (r) PetscValidHeaderSpecific(r,VEC_CLASSID,2); 657 658 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 659 ierr = DMTSSetRHSFunction(dm,f,ctx);CHKERRQ(ierr); 660 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 661 if (!r && !ts->dm && ts->vec_sol) { 662 ierr = VecDuplicate(ts->vec_sol,&ralloc);CHKERRQ(ierr); 663 r = ralloc; 664 } 665 ierr = SNESSetFunction(snes,r,SNESTSFormFunction,ts);CHKERRQ(ierr); 666 ierr = VecDestroy(&ralloc);CHKERRQ(ierr); 667 PetscFunctionReturn(0); 668 } 669 670 #undef __FUNCT__ 671 #define __FUNCT__ "TSSetSolutionFunction" 672 /*@C 673 TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE 674 675 Logically Collective on TS 676 677 Input Parameters: 678 + ts - the TS context obtained from TSCreate() 679 . f - routine for evaluating the solution 680 - ctx - [optional] user-defined context for private data for the 681 function evaluation routine (may be PETSC_NULL) 682 683 Calling sequence of func: 684 $ func (TS ts,PetscReal t,Vec u,void *ctx); 685 686 + t - current timestep 687 . u - output vector 688 - ctx - [optional] user-defined function context 689 690 Notes: 691 This routine is used for testing accuracy of time integration schemes when you already know the solution. 692 If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to 693 create closed-form solutions with non-physical forcing terms. 694 695 For low-dimensional problems solved in serial, such as small discrete systems, TSMonitorErrorODE() can be used to monitor the error history. 696 697 Level: beginner 698 699 .keywords: TS, timestep, set, right-hand-side, function 700 701 .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction() 702 @*/ 703 PetscErrorCode TSSetSolutionFunction(TS ts,PetscErrorCode (*f)(TS,PetscReal,Vec,void*),void *ctx) 704 { 705 PetscErrorCode ierr; 706 DM dm; 707 708 PetscFunctionBegin; 709 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 710 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 711 ierr = DMTSSetSolutionFunction(dm,f,ctx);CHKERRQ(ierr); 712 PetscFunctionReturn(0); 713 } 714 715 #undef __FUNCT__ 716 #define __FUNCT__ "TSSetRHSJacobian" 717 /*@C 718 TSSetRHSJacobian - Sets the function to compute the Jacobian of F, 719 where U_t = F(U,t), as well as the location to store the matrix. 720 721 Logically Collective on TS 722 723 Input Parameters: 724 + ts - the TS context obtained from TSCreate() 725 . A - Jacobian matrix 726 . B - preconditioner matrix (usually same as A) 727 . f - the Jacobian evaluation routine 728 - ctx - [optional] user-defined context for private data for the 729 Jacobian evaluation routine (may be PETSC_NULL) 730 731 Calling sequence of func: 732 $ func (TS ts,PetscReal t,Vec u,Mat *A,Mat *B,MatStructure *flag,void *ctx); 733 734 + t - current timestep 735 . u - input vector 736 . A - matrix A, where U_t = A(t)u 737 . B - preconditioner matrix, usually the same as A 738 . flag - flag indicating information about the preconditioner matrix 739 structure (same as flag in KSPSetOperators()) 740 - ctx - [optional] user-defined context for matrix evaluation routine 741 742 Notes: 743 See KSPSetOperators() for important information about setting the flag 744 output parameter in the routine func(). Be sure to read this information! 745 746 The routine func() takes Mat * as the matrix arguments rather than Mat. 747 This allows the matrix evaluation routine to replace A and/or B with a 748 completely new matrix structure (not just different matrix elements) 749 when appropriate, for instance, if the nonzero structure is changing 750 throughout the global iterations. 751 752 Level: beginner 753 754 .keywords: TS, timestep, set, right-hand-side, Jacobian 755 756 .seealso: SNESDefaultComputeJacobianColor(), TSSetRHSFunction() 757 758 @*/ 759 PetscErrorCode TSSetRHSJacobian(TS ts,Mat A,Mat B,TSRHSJacobian f,void *ctx) 760 { 761 PetscErrorCode ierr; 762 SNES snes; 763 DM dm; 764 TSIJacobian ijacobian; 765 766 PetscFunctionBegin; 767 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 768 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,2); 769 if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3); 770 if (A) PetscCheckSameComm(ts,1,A,2); 771 if (B) PetscCheckSameComm(ts,1,B,3); 772 773 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 774 ierr = DMTSSetRHSJacobian(dm,f,ctx);CHKERRQ(ierr); 775 ierr = DMTSGetIJacobian(dm,&ijacobian,PETSC_NULL);CHKERRQ(ierr); 776 777 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 778 if (!ijacobian) { 779 ierr = SNESSetJacobian(snes,A,B,SNESTSFormJacobian,ts);CHKERRQ(ierr); 780 } 781 if (A) { 782 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 783 ierr = MatDestroy(&ts->Arhs);CHKERRQ(ierr); 784 ts->Arhs = A; 785 } 786 if (B) { 787 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr); 788 ierr = MatDestroy(&ts->Brhs);CHKERRQ(ierr); 789 ts->Brhs = B; 790 } 791 PetscFunctionReturn(0); 792 } 793 794 795 #undef __FUNCT__ 796 #define __FUNCT__ "TSSetIFunction" 797 /*@C 798 TSSetIFunction - Set the function to compute F(t,U,U_t) where F = 0 is the DAE to be solved. 799 800 Logically Collective on TS 801 802 Input Parameters: 803 + ts - the TS context obtained from TSCreate() 804 . r - vector to hold the residual (or PETSC_NULL to have it created internally) 805 . f - the function evaluation routine 806 - ctx - user-defined context for private data for the function evaluation routine (may be PETSC_NULL) 807 808 Calling sequence of f: 809 $ f(TS ts,PetscReal t,Vec u,Vec u_t,Vec F,ctx); 810 811 + t - time at step/stage being solved 812 . u - state vector 813 . u_t - time derivative of state vector 814 . F - function vector 815 - ctx - [optional] user-defined context for matrix evaluation routine 816 817 Important: 818 The user MUST call either this routine, TSSetRHSFunction(). This routine must be used when not solving an ODE, for example a DAE. 819 820 Level: beginner 821 822 .keywords: TS, timestep, set, DAE, Jacobian 823 824 .seealso: TSSetRHSJacobian(), TSSetRHSFunction(), TSSetIJacobian() 825 @*/ 826 PetscErrorCode TSSetIFunction(TS ts,Vec res,TSIFunction f,void *ctx) 827 { 828 PetscErrorCode ierr; 829 SNES snes; 830 Vec resalloc = PETSC_NULL; 831 DM dm; 832 833 PetscFunctionBegin; 834 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 835 if (res) PetscValidHeaderSpecific(res,VEC_CLASSID,2); 836 837 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 838 ierr = DMTSSetIFunction(dm,f,ctx);CHKERRQ(ierr); 839 840 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 841 if (!res && !ts->dm && ts->vec_sol) { 842 ierr = VecDuplicate(ts->vec_sol,&resalloc);CHKERRQ(ierr); 843 res = resalloc; 844 } 845 ierr = SNESSetFunction(snes,res,SNESTSFormFunction,ts);CHKERRQ(ierr); 846 ierr = VecDestroy(&resalloc);CHKERRQ(ierr); 847 848 PetscFunctionReturn(0); 849 } 850 851 #undef __FUNCT__ 852 #define __FUNCT__ "TSGetIFunction" 853 /*@C 854 TSGetIFunction - Returns the vector where the implicit residual is stored and the function/contex to compute it. 855 856 Not Collective 857 858 Input Parameter: 859 . ts - the TS context 860 861 Output Parameter: 862 + r - vector to hold residual (or PETSC_NULL) 863 . func - the function to compute residual (or PETSC_NULL) 864 - ctx - the function context (or PETSC_NULL) 865 866 Level: advanced 867 868 .keywords: TS, nonlinear, get, function 869 870 .seealso: TSSetIFunction(), SNESGetFunction() 871 @*/ 872 PetscErrorCode TSGetIFunction(TS ts,Vec *r,TSIFunction *func,void **ctx) 873 { 874 PetscErrorCode ierr; 875 SNES snes; 876 DM dm; 877 878 PetscFunctionBegin; 879 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 880 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 881 ierr = SNESGetFunction(snes,r,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 882 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 883 ierr = DMTSGetIFunction(dm,func,ctx);CHKERRQ(ierr); 884 PetscFunctionReturn(0); 885 } 886 887 #undef __FUNCT__ 888 #define __FUNCT__ "TSGetRHSFunction" 889 /*@C 890 TSGetRHSFunction - Returns the vector where the right hand side is stored and the function/context to compute it. 891 892 Not Collective 893 894 Input Parameter: 895 . ts - the TS context 896 897 Output Parameter: 898 + r - vector to hold computed right hand side (or PETSC_NULL) 899 . func - the function to compute right hand side (or PETSC_NULL) 900 - ctx - the function context (or PETSC_NULL) 901 902 Level: advanced 903 904 .keywords: TS, nonlinear, get, function 905 906 .seealso: TSSetRhsfunction(), SNESGetFunction() 907 @*/ 908 PetscErrorCode TSGetRHSFunction(TS ts,Vec *r,TSRHSFunction *func,void **ctx) 909 { 910 PetscErrorCode ierr; 911 SNES snes; 912 DM dm; 913 914 PetscFunctionBegin; 915 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 916 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 917 ierr = SNESGetFunction(snes,r,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 918 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 919 ierr = DMTSGetRHSFunction(dm,func,ctx);CHKERRQ(ierr); 920 PetscFunctionReturn(0); 921 } 922 923 #undef __FUNCT__ 924 #define __FUNCT__ "TSSetIJacobian" 925 /*@C 926 TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function 927 you provided with TSSetIFunction(). 928 929 Logically Collective on TS 930 931 Input Parameters: 932 + ts - the TS context obtained from TSCreate() 933 . A - Jacobian matrix 934 . B - preconditioning matrix for A (may be same as A) 935 . f - the Jacobian evaluation routine 936 - ctx - user-defined context for private data for the Jacobian evaluation routine (may be PETSC_NULL) 937 938 Calling sequence of f: 939 $ f(TS ts,PetscReal t,Vec U,Vec U_t,PetscReal a,Mat *A,Mat *B,MatStructure *flag,void *ctx); 940 941 + t - time at step/stage being solved 942 . U - state vector 943 . U_t - time derivative of state vector 944 . a - shift 945 . A - Jacobian of G(U) = F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t 946 . B - preconditioning matrix for A, may be same as A 947 . flag - flag indicating information about the preconditioner matrix 948 structure (same as flag in KSPSetOperators()) 949 - ctx - [optional] user-defined context for matrix evaluation routine 950 951 Notes: 952 The matrices A and B are exactly the matrices that are used by SNES for the nonlinear solve. 953 954 The matrix dF/dU + a*dF/dU_t you provide turns out to be 955 the Jacobian of G(U) = F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved. 956 The time integrator internally approximates U_t by W+a*U where the positive "shift" 957 a and vector W depend on the integration method, step size, and past states. For example with 958 the backward Euler method a = 1/dt and W = -a*U(previous timestep) so 959 W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt 960 961 Level: beginner 962 963 .keywords: TS, timestep, DAE, Jacobian 964 965 .seealso: TSSetIFunction(), TSSetRHSJacobian(), SNESDefaultComputeJacobianColor(), SNESDefaultComputeJacobian() 966 967 @*/ 968 PetscErrorCode TSSetIJacobian(TS ts,Mat A,Mat B,TSIJacobian f,void *ctx) 969 { 970 PetscErrorCode ierr; 971 SNES snes; 972 DM dm; 973 974 PetscFunctionBegin; 975 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 976 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,2); 977 if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3); 978 if (A) PetscCheckSameComm(ts,1,A,2); 979 if (B) PetscCheckSameComm(ts,1,B,3); 980 981 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 982 ierr = DMTSSetIJacobian(dm,f,ctx);CHKERRQ(ierr); 983 984 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 985 ierr = SNESSetJacobian(snes,A,B,SNESTSFormJacobian,ts);CHKERRQ(ierr); 986 PetscFunctionReturn(0); 987 } 988 989 #undef __FUNCT__ 990 #define __FUNCT__ "TSView" 991 /*@C 992 TSView - Prints the TS data structure. 993 994 Collective on TS 995 996 Input Parameters: 997 + ts - the TS context obtained from TSCreate() 998 - viewer - visualization context 999 1000 Options Database Key: 1001 . -ts_view - calls TSView() at end of TSStep() 1002 1003 Notes: 1004 The available visualization contexts include 1005 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 1006 - PETSC_VIEWER_STDOUT_WORLD - synchronized standard 1007 output where only the first processor opens 1008 the file. All other processors send their 1009 data to the first processor to print. 1010 1011 The user can open an alternative visualization context with 1012 PetscViewerASCIIOpen() - output to a specified file. 1013 1014 Level: beginner 1015 1016 .keywords: TS, timestep, view 1017 1018 .seealso: PetscViewerASCIIOpen() 1019 @*/ 1020 PetscErrorCode TSView(TS ts,PetscViewer viewer) 1021 { 1022 PetscErrorCode ierr; 1023 TSType type; 1024 PetscBool iascii,isstring,isundials; 1025 1026 PetscFunctionBegin; 1027 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1028 if (!viewer) { 1029 ierr = PetscViewerASCIIGetStdout(((PetscObject)ts)->comm,&viewer);CHKERRQ(ierr); 1030 } 1031 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 1032 PetscCheckSameComm(ts,1,viewer,2); 1033 1034 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1035 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 1036 if (iascii) { 1037 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)ts,viewer,"TS Object");CHKERRQ(ierr); 1038 ierr = PetscViewerASCIIPrintf(viewer," maximum steps=%D\n",ts->max_steps);CHKERRQ(ierr); 1039 ierr = PetscViewerASCIIPrintf(viewer," maximum time=%G\n",ts->max_time);CHKERRQ(ierr); 1040 if (ts->problem_type == TS_NONLINEAR) { 1041 ierr = PetscViewerASCIIPrintf(viewer," total number of nonlinear solver iterations=%D\n",ts->snes_its);CHKERRQ(ierr); 1042 ierr = PetscViewerASCIIPrintf(viewer," total number of nonlinear solve failures=%D\n",ts->num_snes_failures);CHKERRQ(ierr); 1043 } 1044 ierr = PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",ts->ksp_its);CHKERRQ(ierr); 1045 ierr = PetscViewerASCIIPrintf(viewer," total number of rejected steps=%D\n",ts->reject);CHKERRQ(ierr); 1046 if (ts->ops->view) { 1047 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1048 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1049 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1050 } 1051 } else if (isstring) { 1052 ierr = TSGetType(ts,&type);CHKERRQ(ierr); 1053 ierr = PetscViewerStringSPrintf(viewer," %-7.7s",type);CHKERRQ(ierr); 1054 } 1055 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1056 ierr = PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&isundials);CHKERRQ(ierr); 1057 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1058 PetscFunctionReturn(0); 1059 } 1060 1061 1062 #undef __FUNCT__ 1063 #define __FUNCT__ "TSSetApplicationContext" 1064 /*@ 1065 TSSetApplicationContext - Sets an optional user-defined context for 1066 the timesteppers. 1067 1068 Logically Collective on TS 1069 1070 Input Parameters: 1071 + ts - the TS context obtained from TSCreate() 1072 - usrP - optional user context 1073 1074 Level: intermediate 1075 1076 .keywords: TS, timestep, set, application, context 1077 1078 .seealso: TSGetApplicationContext() 1079 @*/ 1080 PetscErrorCode TSSetApplicationContext(TS ts,void *usrP) 1081 { 1082 PetscFunctionBegin; 1083 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1084 ts->user = usrP; 1085 PetscFunctionReturn(0); 1086 } 1087 1088 #undef __FUNCT__ 1089 #define __FUNCT__ "TSGetApplicationContext" 1090 /*@ 1091 TSGetApplicationContext - Gets the user-defined context for the 1092 timestepper. 1093 1094 Not Collective 1095 1096 Input Parameter: 1097 . ts - the TS context obtained from TSCreate() 1098 1099 Output Parameter: 1100 . usrP - user context 1101 1102 Level: intermediate 1103 1104 .keywords: TS, timestep, get, application, context 1105 1106 .seealso: TSSetApplicationContext() 1107 @*/ 1108 PetscErrorCode TSGetApplicationContext(TS ts,void *usrP) 1109 { 1110 PetscFunctionBegin; 1111 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1112 *(void**)usrP = ts->user; 1113 PetscFunctionReturn(0); 1114 } 1115 1116 #undef __FUNCT__ 1117 #define __FUNCT__ "TSGetTimeStepNumber" 1118 /*@ 1119 TSGetTimeStepNumber - Gets the number of time steps completed. 1120 1121 Not Collective 1122 1123 Input Parameter: 1124 . ts - the TS context obtained from TSCreate() 1125 1126 Output Parameter: 1127 . iter - number of steps completed so far 1128 1129 Level: intermediate 1130 1131 .keywords: TS, timestep, get, iteration, number 1132 .seealso: TSGetTime(), TSGetTimeStep(), TSSetPreStep(), TSSetPreStage(), TSSetPostStep() 1133 @*/ 1134 PetscErrorCode TSGetTimeStepNumber(TS ts,PetscInt* iter) 1135 { 1136 PetscFunctionBegin; 1137 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1138 PetscValidIntPointer(iter,2); 1139 *iter = ts->steps; 1140 PetscFunctionReturn(0); 1141 } 1142 1143 #undef __FUNCT__ 1144 #define __FUNCT__ "TSSetInitialTimeStep" 1145 /*@ 1146 TSSetInitialTimeStep - Sets the initial timestep to be used, 1147 as well as the initial time. 1148 1149 Logically Collective on TS 1150 1151 Input Parameters: 1152 + ts - the TS context obtained from TSCreate() 1153 . initial_time - the initial time 1154 - time_step - the size of the timestep 1155 1156 Level: intermediate 1157 1158 .seealso: TSSetTimeStep(), TSGetTimeStep() 1159 1160 .keywords: TS, set, initial, timestep 1161 @*/ 1162 PetscErrorCode TSSetInitialTimeStep(TS ts,PetscReal initial_time,PetscReal time_step) 1163 { 1164 PetscErrorCode ierr; 1165 1166 PetscFunctionBegin; 1167 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1168 ierr = TSSetTimeStep(ts,time_step);CHKERRQ(ierr); 1169 ierr = TSSetTime(ts,initial_time);CHKERRQ(ierr); 1170 PetscFunctionReturn(0); 1171 } 1172 1173 #undef __FUNCT__ 1174 #define __FUNCT__ "TSSetTimeStep" 1175 /*@ 1176 TSSetTimeStep - Allows one to reset the timestep at any time, 1177 useful for simple pseudo-timestepping codes. 1178 1179 Logically Collective on TS 1180 1181 Input Parameters: 1182 + ts - the TS context obtained from TSCreate() 1183 - time_step - the size of the timestep 1184 1185 Level: intermediate 1186 1187 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 1188 1189 .keywords: TS, set, timestep 1190 @*/ 1191 PetscErrorCode TSSetTimeStep(TS ts,PetscReal time_step) 1192 { 1193 PetscFunctionBegin; 1194 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1195 PetscValidLogicalCollectiveReal(ts,time_step,2); 1196 ts->time_step = time_step; 1197 ts->time_step_orig = time_step; 1198 PetscFunctionReturn(0); 1199 } 1200 1201 #undef __FUNCT__ 1202 #define __FUNCT__ "TSSetExactFinalTime" 1203 /*@ 1204 TSSetExactFinalTime - Determines whether to interpolate solution to the 1205 exact final time requested by the user or just returns it at the final time 1206 it computed. 1207 1208 Logically Collective on TS 1209 1210 Input Parameter: 1211 + ts - the time-step context 1212 - ft - PETSC_TRUE if interpolates, else PETSC_FALSE 1213 1214 Level: beginner 1215 1216 .seealso: TSSetDuration() 1217 @*/ 1218 PetscErrorCode TSSetExactFinalTime(TS ts,PetscBool flg) 1219 { 1220 1221 PetscFunctionBegin; 1222 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1223 PetscValidLogicalCollectiveBool(ts,flg,2); 1224 ts->exact_final_time = flg; 1225 PetscFunctionReturn(0); 1226 } 1227 1228 #undef __FUNCT__ 1229 #define __FUNCT__ "TSGetTimeStep" 1230 /*@ 1231 TSGetTimeStep - Gets the current timestep size. 1232 1233 Not Collective 1234 1235 Input Parameter: 1236 . ts - the TS context obtained from TSCreate() 1237 1238 Output Parameter: 1239 . dt - the current timestep size 1240 1241 Level: intermediate 1242 1243 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 1244 1245 .keywords: TS, get, timestep 1246 @*/ 1247 PetscErrorCode TSGetTimeStep(TS ts,PetscReal* dt) 1248 { 1249 PetscFunctionBegin; 1250 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1251 PetscValidRealPointer(dt,2); 1252 *dt = ts->time_step; 1253 PetscFunctionReturn(0); 1254 } 1255 1256 #undef __FUNCT__ 1257 #define __FUNCT__ "TSGetSolution" 1258 /*@ 1259 TSGetSolution - Returns the solution at the present timestep. It 1260 is valid to call this routine inside the function that you are evaluating 1261 in order to move to the new timestep. This vector not changed until 1262 the solution at the next timestep has been calculated. 1263 1264 Not Collective, but Vec returned is parallel if TS is parallel 1265 1266 Input Parameter: 1267 . ts - the TS context obtained from TSCreate() 1268 1269 Output Parameter: 1270 . v - the vector containing the solution 1271 1272 Level: intermediate 1273 1274 .seealso: TSGetTimeStep() 1275 1276 .keywords: TS, timestep, get, solution 1277 @*/ 1278 PetscErrorCode TSGetSolution(TS ts,Vec *v) 1279 { 1280 PetscFunctionBegin; 1281 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1282 PetscValidPointer(v,2); 1283 *v = ts->vec_sol; 1284 PetscFunctionReturn(0); 1285 } 1286 1287 /* ----- Routines to initialize and destroy a timestepper ---- */ 1288 #undef __FUNCT__ 1289 #define __FUNCT__ "TSSetProblemType" 1290 /*@ 1291 TSSetProblemType - Sets the type of problem to be solved. 1292 1293 Not collective 1294 1295 Input Parameters: 1296 + ts - The TS 1297 - type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms 1298 .vb 1299 U_t = A U 1300 U_t = A(t) U 1301 U_t = F(t,U) 1302 .ve 1303 1304 Level: beginner 1305 1306 .keywords: TS, problem type 1307 .seealso: TSSetUp(), TSProblemType, TS 1308 @*/ 1309 PetscErrorCode TSSetProblemType(TS ts, TSProblemType type) 1310 { 1311 PetscErrorCode ierr; 1312 1313 PetscFunctionBegin; 1314 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1315 ts->problem_type = type; 1316 if (type == TS_LINEAR) { 1317 SNES snes; 1318 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1319 ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 1320 } 1321 PetscFunctionReturn(0); 1322 } 1323 1324 #undef __FUNCT__ 1325 #define __FUNCT__ "TSGetProblemType" 1326 /*@C 1327 TSGetProblemType - Gets the type of problem to be solved. 1328 1329 Not collective 1330 1331 Input Parameter: 1332 . ts - The TS 1333 1334 Output Parameter: 1335 . type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms 1336 .vb 1337 M U_t = A U 1338 M(t) U_t = A(t) U 1339 U_t = F(t,U) 1340 .ve 1341 1342 Level: beginner 1343 1344 .keywords: TS, problem type 1345 .seealso: TSSetUp(), TSProblemType, TS 1346 @*/ 1347 PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type) 1348 { 1349 PetscFunctionBegin; 1350 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1351 PetscValidIntPointer(type,2); 1352 *type = ts->problem_type; 1353 PetscFunctionReturn(0); 1354 } 1355 1356 #undef __FUNCT__ 1357 #define __FUNCT__ "TSSetUp" 1358 /*@ 1359 TSSetUp - Sets up the internal data structures for the later use 1360 of a timestepper. 1361 1362 Collective on TS 1363 1364 Input Parameter: 1365 . ts - the TS context obtained from TSCreate() 1366 1367 Notes: 1368 For basic use of the TS solvers the user need not explicitly call 1369 TSSetUp(), since these actions will automatically occur during 1370 the call to TSStep(). However, if one wishes to control this 1371 phase separately, TSSetUp() should be called after TSCreate() 1372 and optional routines of the form TSSetXXX(), but before TSStep(). 1373 1374 Level: advanced 1375 1376 .keywords: TS, timestep, setup 1377 1378 .seealso: TSCreate(), TSStep(), TSDestroy() 1379 @*/ 1380 PetscErrorCode TSSetUp(TS ts) 1381 { 1382 PetscErrorCode ierr; 1383 DM dm; 1384 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 1385 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 1386 TSIJacobian ijac; 1387 TSRHSJacobian rhsjac; 1388 1389 PetscFunctionBegin; 1390 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1391 if (ts->setupcalled) PetscFunctionReturn(0); 1392 1393 if (!((PetscObject)ts)->type_name) { 1394 ierr = TSSetType(ts,TSEULER);CHKERRQ(ierr); 1395 } 1396 if (ts->exact_final_time == PETSC_DECIDE) ts->exact_final_time = PETSC_FALSE; 1397 1398 if (!ts->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetSolution() first"); 1399 1400 ierr = TSGetAdapt(ts,&ts->adapt);CHKERRQ(ierr); 1401 1402 if (ts->ops->setup) { 1403 ierr = (*ts->ops->setup)(ts);CHKERRQ(ierr); 1404 } 1405 1406 /* in the case where we've set a DMTSFunction or what have you, we need the default SNESFunction 1407 to be set right but can't do it elsewhere due to the overreliance on ctx=ts. 1408 */ 1409 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1410 ierr = DMSNESGetFunction(dm,&func,PETSC_NULL);CHKERRQ(ierr); 1411 if (!func) { 1412 ierr =DMSNESSetFunction(dm,SNESTSFormFunction,ts);CHKERRQ(ierr); 1413 } 1414 /* if the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it. 1415 Otherwise, the SNES will use coloring internally to form the Jacobian. 1416 */ 1417 ierr = DMSNESGetJacobian(dm,&jac,PETSC_NULL);CHKERRQ(ierr); 1418 ierr = DMTSGetIJacobian(dm,&ijac,PETSC_NULL);CHKERRQ(ierr); 1419 ierr = DMTSGetRHSJacobian(dm,&rhsjac,PETSC_NULL);CHKERRQ(ierr); 1420 if (!jac && (ijac || rhsjac)) { 1421 ierr = DMSNESSetJacobian(dm,SNESTSFormJacobian,ts);CHKERRQ(ierr); 1422 } 1423 ts->setupcalled = PETSC_TRUE; 1424 PetscFunctionReturn(0); 1425 } 1426 1427 #undef __FUNCT__ 1428 #define __FUNCT__ "TSReset" 1429 /*@ 1430 TSReset - Resets a TS context and removes any allocated Vecs and Mats. 1431 1432 Collective on TS 1433 1434 Input Parameter: 1435 . ts - the TS context obtained from TSCreate() 1436 1437 Level: beginner 1438 1439 .keywords: TS, timestep, reset 1440 1441 .seealso: TSCreate(), TSSetup(), TSDestroy() 1442 @*/ 1443 PetscErrorCode TSReset(TS ts) 1444 { 1445 PetscErrorCode ierr; 1446 1447 PetscFunctionBegin; 1448 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1449 if (ts->ops->reset) { 1450 ierr = (*ts->ops->reset)(ts);CHKERRQ(ierr); 1451 } 1452 if (ts->snes) {ierr = SNESReset(ts->snes);CHKERRQ(ierr);} 1453 ierr = MatDestroy(&ts->Arhs);CHKERRQ(ierr); 1454 ierr = MatDestroy(&ts->Brhs);CHKERRQ(ierr); 1455 ierr = VecDestroy(&ts->Frhs);CHKERRQ(ierr); 1456 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 1457 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 1458 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 1459 ierr = VecDestroyVecs(ts->nwork,&ts->work);CHKERRQ(ierr); 1460 ts->setupcalled = PETSC_FALSE; 1461 PetscFunctionReturn(0); 1462 } 1463 1464 #undef __FUNCT__ 1465 #define __FUNCT__ "TSDestroy" 1466 /*@ 1467 TSDestroy - Destroys the timestepper context that was created 1468 with TSCreate(). 1469 1470 Collective on TS 1471 1472 Input Parameter: 1473 . ts - the TS context obtained from TSCreate() 1474 1475 Level: beginner 1476 1477 .keywords: TS, timestepper, destroy 1478 1479 .seealso: TSCreate(), TSSetUp(), TSSolve() 1480 @*/ 1481 PetscErrorCode TSDestroy(TS *ts) 1482 { 1483 PetscErrorCode ierr; 1484 1485 PetscFunctionBegin; 1486 if (!*ts) PetscFunctionReturn(0); 1487 PetscValidHeaderSpecific((*ts),TS_CLASSID,1); 1488 if (--((PetscObject)(*ts))->refct > 0) {*ts = 0; PetscFunctionReturn(0);} 1489 1490 ierr = TSReset((*ts));CHKERRQ(ierr); 1491 1492 /* if memory was published with AMS then destroy it */ 1493 ierr = PetscObjectDepublish((*ts));CHKERRQ(ierr); 1494 if ((*ts)->ops->destroy) {ierr = (*(*ts)->ops->destroy)((*ts));CHKERRQ(ierr);} 1495 1496 ierr = TSAdaptDestroy(&(*ts)->adapt);CHKERRQ(ierr); 1497 ierr = SNESDestroy(&(*ts)->snes);CHKERRQ(ierr); 1498 ierr = DMDestroy(&(*ts)->dm);CHKERRQ(ierr); 1499 ierr = TSMonitorCancel((*ts));CHKERRQ(ierr); 1500 1501 ierr = PetscHeaderDestroy(ts);CHKERRQ(ierr); 1502 PetscFunctionReturn(0); 1503 } 1504 1505 #undef __FUNCT__ 1506 #define __FUNCT__ "TSGetSNES" 1507 /*@ 1508 TSGetSNES - Returns the SNES (nonlinear solver) associated with 1509 a TS (timestepper) context. Valid only for nonlinear problems. 1510 1511 Not Collective, but SNES is parallel if TS is parallel 1512 1513 Input Parameter: 1514 . ts - the TS context obtained from TSCreate() 1515 1516 Output Parameter: 1517 . snes - the nonlinear solver context 1518 1519 Notes: 1520 The user can then directly manipulate the SNES context to set various 1521 options, etc. Likewise, the user can then extract and manipulate the 1522 KSP, KSP, and PC contexts as well. 1523 1524 TSGetSNES() does not work for integrators that do not use SNES; in 1525 this case TSGetSNES() returns PETSC_NULL in snes. 1526 1527 Level: beginner 1528 1529 .keywords: timestep, get, SNES 1530 @*/ 1531 PetscErrorCode TSGetSNES(TS ts,SNES *snes) 1532 { 1533 PetscErrorCode ierr; 1534 1535 PetscFunctionBegin; 1536 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1537 PetscValidPointer(snes,2); 1538 if (!ts->snes) { 1539 ierr = SNESCreate(((PetscObject)ts)->comm,&ts->snes);CHKERRQ(ierr); 1540 ierr = SNESSetFunction(ts->snes,PETSC_NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 1541 ierr = PetscLogObjectParent(ts,ts->snes);CHKERRQ(ierr); 1542 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);CHKERRQ(ierr); 1543 if (ts->dm) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 1544 if (ts->problem_type == TS_LINEAR) { 1545 ierr = SNESSetType(ts->snes,SNESKSPONLY);CHKERRQ(ierr); 1546 } 1547 } 1548 *snes = ts->snes; 1549 PetscFunctionReturn(0); 1550 } 1551 1552 #undef __FUNCT__ 1553 #define __FUNCT__ "TSGetKSP" 1554 /*@ 1555 TSGetKSP - Returns the KSP (linear solver) associated with 1556 a TS (timestepper) context. 1557 1558 Not Collective, but KSP is parallel if TS is parallel 1559 1560 Input Parameter: 1561 . ts - the TS context obtained from TSCreate() 1562 1563 Output Parameter: 1564 . ksp - the nonlinear solver context 1565 1566 Notes: 1567 The user can then directly manipulate the KSP context to set various 1568 options, etc. Likewise, the user can then extract and manipulate the 1569 KSP and PC contexts as well. 1570 1571 TSGetKSP() does not work for integrators that do not use KSP; 1572 in this case TSGetKSP() returns PETSC_NULL in ksp. 1573 1574 Level: beginner 1575 1576 .keywords: timestep, get, KSP 1577 @*/ 1578 PetscErrorCode TSGetKSP(TS ts,KSP *ksp) 1579 { 1580 PetscErrorCode ierr; 1581 SNES snes; 1582 1583 PetscFunctionBegin; 1584 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1585 PetscValidPointer(ksp,2); 1586 if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first"); 1587 if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()"); 1588 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1589 ierr = SNESGetKSP(snes,ksp);CHKERRQ(ierr); 1590 PetscFunctionReturn(0); 1591 } 1592 1593 /* ----------- Routines to set solver parameters ---------- */ 1594 1595 #undef __FUNCT__ 1596 #define __FUNCT__ "TSGetDuration" 1597 /*@ 1598 TSGetDuration - Gets the maximum number of timesteps to use and 1599 maximum time for iteration. 1600 1601 Not Collective 1602 1603 Input Parameters: 1604 + ts - the TS context obtained from TSCreate() 1605 . maxsteps - maximum number of iterations to use, or PETSC_NULL 1606 - maxtime - final time to iterate to, or PETSC_NULL 1607 1608 Level: intermediate 1609 1610 .keywords: TS, timestep, get, maximum, iterations, time 1611 @*/ 1612 PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime) 1613 { 1614 PetscFunctionBegin; 1615 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1616 if (maxsteps) { 1617 PetscValidIntPointer(maxsteps,2); 1618 *maxsteps = ts->max_steps; 1619 } 1620 if (maxtime) { 1621 PetscValidScalarPointer(maxtime,3); 1622 *maxtime = ts->max_time; 1623 } 1624 PetscFunctionReturn(0); 1625 } 1626 1627 #undef __FUNCT__ 1628 #define __FUNCT__ "TSSetDuration" 1629 /*@ 1630 TSSetDuration - Sets the maximum number of timesteps to use and 1631 maximum time for iteration. 1632 1633 Logically Collective on TS 1634 1635 Input Parameters: 1636 + ts - the TS context obtained from TSCreate() 1637 . maxsteps - maximum number of iterations to use 1638 - maxtime - final time to iterate to 1639 1640 Options Database Keys: 1641 . -ts_max_steps <maxsteps> - Sets maxsteps 1642 . -ts_final_time <maxtime> - Sets maxtime 1643 1644 Notes: 1645 The default maximum number of iterations is 5000. Default time is 5.0 1646 1647 Level: intermediate 1648 1649 .keywords: TS, timestep, set, maximum, iterations 1650 1651 .seealso: TSSetExactFinalTime() 1652 @*/ 1653 PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime) 1654 { 1655 PetscFunctionBegin; 1656 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1657 PetscValidLogicalCollectiveInt(ts,maxsteps,2); 1658 PetscValidLogicalCollectiveReal(ts,maxtime,2); 1659 if (maxsteps >= 0) ts->max_steps = maxsteps; 1660 if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime; 1661 PetscFunctionReturn(0); 1662 } 1663 1664 #undef __FUNCT__ 1665 #define __FUNCT__ "TSSetSolution" 1666 /*@ 1667 TSSetSolution - Sets the initial solution vector 1668 for use by the TS routines. 1669 1670 Logically Collective on TS and Vec 1671 1672 Input Parameters: 1673 + ts - the TS context obtained from TSCreate() 1674 - x - the solution vector 1675 1676 Level: beginner 1677 1678 .keywords: TS, timestep, set, solution, initial conditions 1679 @*/ 1680 PetscErrorCode TSSetSolution(TS ts,Vec x) 1681 { 1682 PetscErrorCode ierr; 1683 DM dm; 1684 1685 PetscFunctionBegin; 1686 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1687 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 1688 ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr); 1689 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 1690 ts->vec_sol = x; 1691 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1692 ierr = DMShellSetGlobalVector(dm,x);CHKERRQ(ierr); 1693 PetscFunctionReturn(0); 1694 } 1695 1696 #undef __FUNCT__ 1697 #define __FUNCT__ "TSSetPreStep" 1698 /*@C 1699 TSSetPreStep - Sets the general-purpose function 1700 called once at the beginning of each time step. 1701 1702 Logically Collective on TS 1703 1704 Input Parameters: 1705 + ts - The TS context obtained from TSCreate() 1706 - func - The function 1707 1708 Calling sequence of func: 1709 . func (TS ts); 1710 1711 Level: intermediate 1712 1713 Note: 1714 If a step is rejected, TSStep() will call this routine again before each attempt. 1715 The last completed time step number can be queried using TSGetTimeStepNumber(), the 1716 size of the step being attempted can be obtained using TSGetTimeStep(). 1717 1718 .keywords: TS, timestep 1719 .seealso: TSSetPreStage(), TSSetPostStep(), TSStep() 1720 @*/ 1721 PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS)) 1722 { 1723 PetscFunctionBegin; 1724 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1725 ts->ops->prestep = func; 1726 PetscFunctionReturn(0); 1727 } 1728 1729 #undef __FUNCT__ 1730 #define __FUNCT__ "TSPreStep" 1731 /*@ 1732 TSPreStep - Runs the user-defined pre-step function. 1733 1734 Collective on TS 1735 1736 Input Parameters: 1737 . ts - The TS context obtained from TSCreate() 1738 1739 Notes: 1740 TSPreStep() is typically used within time stepping implementations, 1741 so most users would not generally call this routine themselves. 1742 1743 Level: developer 1744 1745 .keywords: TS, timestep 1746 .seealso: TSSetPreStep(), TSPreStage(), TSPostStep() 1747 @*/ 1748 PetscErrorCode TSPreStep(TS ts) 1749 { 1750 PetscErrorCode ierr; 1751 1752 PetscFunctionBegin; 1753 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1754 if (ts->ops->prestep) { 1755 PetscStackPush("TS PreStep function"); 1756 ierr = (*ts->ops->prestep)(ts);CHKERRQ(ierr); 1757 PetscStackPop; 1758 } 1759 PetscFunctionReturn(0); 1760 } 1761 1762 #undef __FUNCT__ 1763 #define __FUNCT__ "TSSetPreStage" 1764 /*@C 1765 TSSetPreStage - Sets the general-purpose function 1766 called once at the beginning of each stage. 1767 1768 Logically Collective on TS 1769 1770 Input Parameters: 1771 + ts - The TS context obtained from TSCreate() 1772 - func - The function 1773 1774 Calling sequence of func: 1775 . PetscErrorCode func(TS ts, PetscReal stagetime); 1776 1777 Level: intermediate 1778 1779 Note: 1780 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 1781 The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being 1782 attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime(). 1783 1784 .keywords: TS, timestep 1785 .seealso: TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 1786 @*/ 1787 PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS,PetscReal)) 1788 { 1789 PetscFunctionBegin; 1790 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1791 ts->ops->prestage = func; 1792 PetscFunctionReturn(0); 1793 } 1794 1795 #undef __FUNCT__ 1796 #define __FUNCT__ "TSPreStage" 1797 /*@ 1798 TSPreStage - Runs the user-defined pre-stage function set using TSSetPreStage() 1799 1800 Collective on TS 1801 1802 Input Parameters: 1803 . ts - The TS context obtained from TSCreate() 1804 1805 Notes: 1806 TSPreStage() is typically used within time stepping implementations, 1807 most users would not generally call this routine themselves. 1808 1809 Level: developer 1810 1811 .keywords: TS, timestep 1812 .seealso: TSSetPreStep(), TSPreStep(), TSPostStep() 1813 @*/ 1814 PetscErrorCode TSPreStage(TS ts, PetscReal stagetime) 1815 { 1816 PetscErrorCode ierr; 1817 1818 PetscFunctionBegin; 1819 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1820 if (ts->ops->prestage) { 1821 PetscStackPush("TS PreStage function"); 1822 ierr = (*ts->ops->prestage)(ts,stagetime);CHKERRQ(ierr); 1823 PetscStackPop; 1824 } 1825 PetscFunctionReturn(0); 1826 } 1827 1828 #undef __FUNCT__ 1829 #define __FUNCT__ "TSSetPostStep" 1830 /*@C 1831 TSSetPostStep - Sets the general-purpose function 1832 called once at the end of each time step. 1833 1834 Logically Collective on TS 1835 1836 Input Parameters: 1837 + ts - The TS context obtained from TSCreate() 1838 - func - The function 1839 1840 Calling sequence of func: 1841 $ func (TS ts); 1842 1843 Level: intermediate 1844 1845 .keywords: TS, timestep 1846 .seealso: TSSetPreStep(), TSSetPreStage(), TSGetTimeStep(), TSGetTimeStepNumber(), TSGetTime() 1847 @*/ 1848 PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS)) 1849 { 1850 PetscFunctionBegin; 1851 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1852 ts->ops->poststep = func; 1853 PetscFunctionReturn(0); 1854 } 1855 1856 #undef __FUNCT__ 1857 #define __FUNCT__ "TSPostStep" 1858 /*@ 1859 TSPostStep - Runs the user-defined post-step function. 1860 1861 Collective on TS 1862 1863 Input Parameters: 1864 . ts - The TS context obtained from TSCreate() 1865 1866 Notes: 1867 TSPostStep() is typically used within time stepping implementations, 1868 so most users would not generally call this routine themselves. 1869 1870 Level: developer 1871 1872 .keywords: TS, timestep 1873 @*/ 1874 PetscErrorCode TSPostStep(TS ts) 1875 { 1876 PetscErrorCode ierr; 1877 1878 PetscFunctionBegin; 1879 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1880 if (ts->ops->poststep) { 1881 PetscStackPush("TS PostStep function"); 1882 ierr = (*ts->ops->poststep)(ts);CHKERRQ(ierr); 1883 PetscStackPop; 1884 } 1885 PetscFunctionReturn(0); 1886 } 1887 1888 /* ------------ Routines to set performance monitoring options ----------- */ 1889 1890 #undef __FUNCT__ 1891 #define __FUNCT__ "TSMonitorSet" 1892 /*@C 1893 TSMonitorSet - Sets an ADDITIONAL function that is to be used at every 1894 timestep to display the iteration's progress. 1895 1896 Logically Collective on TS 1897 1898 Input Parameters: 1899 + ts - the TS context obtained from TSCreate() 1900 . monitor - monitoring routine 1901 . mctx - [optional] user-defined context for private data for the 1902 monitor routine (use PETSC_NULL if no context is desired) 1903 - monitordestroy - [optional] routine that frees monitor context 1904 (may be PETSC_NULL) 1905 1906 Calling sequence of monitor: 1907 $ int monitor(TS ts,PetscInt steps,PetscReal time,Vec x,void *mctx) 1908 1909 + ts - the TS context 1910 . steps - iteration number 1911 . time - current time 1912 . x - current iterate 1913 - mctx - [optional] monitoring context 1914 1915 Notes: 1916 This routine adds an additional monitor to the list of monitors that 1917 already has been loaded. 1918 1919 Fortran notes: Only a single monitor function can be set for each TS object 1920 1921 Level: intermediate 1922 1923 .keywords: TS, timestep, set, monitor 1924 1925 .seealso: TSMonitorDefault(), TSMonitorCancel() 1926 @*/ 1927 PetscErrorCode TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**)) 1928 { 1929 PetscFunctionBegin; 1930 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1931 if (ts->numbermonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 1932 ts->monitor[ts->numbermonitors] = monitor; 1933 ts->mdestroy[ts->numbermonitors] = mdestroy; 1934 ts->monitorcontext[ts->numbermonitors++] = (void*)mctx; 1935 PetscFunctionReturn(0); 1936 } 1937 1938 #undef __FUNCT__ 1939 #define __FUNCT__ "TSMonitorCancel" 1940 /*@C 1941 TSMonitorCancel - Clears all the monitors that have been set on a time-step object. 1942 1943 Logically Collective on TS 1944 1945 Input Parameters: 1946 . ts - the TS context obtained from TSCreate() 1947 1948 Notes: 1949 There is no way to remove a single, specific monitor. 1950 1951 Level: intermediate 1952 1953 .keywords: TS, timestep, set, monitor 1954 1955 .seealso: TSMonitorDefault(), TSMonitorSet() 1956 @*/ 1957 PetscErrorCode TSMonitorCancel(TS ts) 1958 { 1959 PetscErrorCode ierr; 1960 PetscInt i; 1961 1962 PetscFunctionBegin; 1963 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1964 for (i=0; i<ts->numbermonitors; i++) { 1965 if (ts->mdestroy[i]) { 1966 ierr = (*ts->mdestroy[i])(&ts->monitorcontext[i]);CHKERRQ(ierr); 1967 } 1968 } 1969 ts->numbermonitors = 0; 1970 PetscFunctionReturn(0); 1971 } 1972 1973 #undef __FUNCT__ 1974 #define __FUNCT__ "TSMonitorDefault" 1975 /*@ 1976 TSMonitorDefault - Sets the Default monitor 1977 1978 Level: intermediate 1979 1980 .keywords: TS, set, monitor 1981 1982 .seealso: TSMonitorDefault(), TSMonitorSet() 1983 @*/ 1984 PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy) 1985 { 1986 PetscErrorCode ierr; 1987 PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(((PetscObject)ts)->comm); 1988 1989 PetscFunctionBegin; 1990 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 1991 ierr = PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g\n",step,(double)ts->time_step,(double)ptime);CHKERRQ(ierr); 1992 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 1993 PetscFunctionReturn(0); 1994 } 1995 1996 #undef __FUNCT__ 1997 #define __FUNCT__ "TSSetRetainStages" 1998 /*@ 1999 TSSetRetainStages - Request that all stages in the upcoming step be stored so that interpolation will be available. 2000 2001 Logically Collective on TS 2002 2003 Input Argument: 2004 . ts - time stepping context 2005 2006 Output Argument: 2007 . flg - PETSC_TRUE or PETSC_FALSE 2008 2009 Level: intermediate 2010 2011 .keywords: TS, set 2012 2013 .seealso: TSInterpolate(), TSSetPostStep() 2014 @*/ 2015 PetscErrorCode TSSetRetainStages(TS ts,PetscBool flg) 2016 { 2017 2018 PetscFunctionBegin; 2019 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2020 ts->retain_stages = flg; 2021 PetscFunctionReturn(0); 2022 } 2023 2024 #undef __FUNCT__ 2025 #define __FUNCT__ "TSInterpolate" 2026 /*@ 2027 TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval 2028 2029 Collective on TS 2030 2031 Input Argument: 2032 + ts - time stepping context 2033 - t - time to interpolate to 2034 2035 Output Argument: 2036 . X - state at given time 2037 2038 Notes: 2039 The user should call TSSetRetainStages() before taking a step in which interpolation will be requested. 2040 2041 Level: intermediate 2042 2043 Developer Notes: 2044 TSInterpolate() and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints. 2045 2046 .keywords: TS, set 2047 2048 .seealso: TSSetRetainStages(), TSSetPostStep() 2049 @*/ 2050 PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec X) 2051 { 2052 PetscErrorCode ierr; 2053 2054 PetscFunctionBegin; 2055 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2056 if (t < ts->ptime - ts->time_step_prev || t > ts->ptime) SETERRQ3(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Requested time %G not in last time steps [%G,%G]",t,ts->ptime-ts->time_step_prev,ts->ptime); 2057 if (!ts->ops->interpolate) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name); 2058 ierr = (*ts->ops->interpolate)(ts,t,X);CHKERRQ(ierr); 2059 PetscFunctionReturn(0); 2060 } 2061 2062 #undef __FUNCT__ 2063 #define __FUNCT__ "TSStep" 2064 /*@ 2065 TSStep - Steps one time step 2066 2067 Collective on TS 2068 2069 Input Parameter: 2070 . ts - the TS context obtained from TSCreate() 2071 2072 Level: intermediate 2073 2074 Notes: 2075 The hook set using TSSetPreStep() is called before each attempt to take the step. In general, the time step size may 2076 be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages. 2077 2078 .keywords: TS, timestep, solve 2079 2080 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage() 2081 @*/ 2082 PetscErrorCode TSStep(TS ts) 2083 { 2084 PetscReal ptime_prev; 2085 PetscErrorCode ierr; 2086 2087 PetscFunctionBegin; 2088 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2089 ierr = TSSetUp(ts);CHKERRQ(ierr); 2090 2091 ts->reason = TS_CONVERGED_ITERATING; 2092 2093 ptime_prev = ts->ptime; 2094 ierr = PetscLogEventBegin(TS_Step,ts,0,0,0);CHKERRQ(ierr); 2095 ierr = (*ts->ops->step)(ts);CHKERRQ(ierr); 2096 ierr = PetscLogEventEnd(TS_Step,ts,0,0,0);CHKERRQ(ierr); 2097 ts->time_step_prev = ts->ptime - ptime_prev; 2098 2099 if (ts->reason < 0) { 2100 if (ts->errorifstepfailed) { 2101 if (ts->reason == TS_DIVERGED_NONLINEAR_SOLVE) { 2102 SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s, increase -ts_max_snes_failures or make negative to attempt recovery",TSConvergedReasons[ts->reason]); 2103 } else SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]); 2104 } 2105 } else if (!ts->reason) { 2106 if (ts->steps >= ts->max_steps) 2107 ts->reason = TS_CONVERGED_ITS; 2108 else if (ts->ptime >= ts->max_time) 2109 ts->reason = TS_CONVERGED_TIME; 2110 } 2111 2112 PetscFunctionReturn(0); 2113 } 2114 2115 #undef __FUNCT__ 2116 #define __FUNCT__ "TSEvaluateStep" 2117 /*@ 2118 TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy. 2119 2120 Collective on TS 2121 2122 Input Arguments: 2123 + ts - time stepping context 2124 . order - desired order of accuracy 2125 - done - whether the step was evaluated at this order (pass PETSC_NULL to generate an error if not available) 2126 2127 Output Arguments: 2128 . X - state at the end of the current step 2129 2130 Level: advanced 2131 2132 Notes: 2133 This function cannot be called until all stages have been evaluated. 2134 It is normally called by adaptive controllers before a step has been accepted and may also be called by the user after TSStep() has returned. 2135 2136 .seealso: TSStep(), TSAdapt 2137 @*/ 2138 PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec X,PetscBool *done) 2139 { 2140 PetscErrorCode ierr; 2141 2142 PetscFunctionBegin; 2143 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2144 PetscValidType(ts,1); 2145 PetscValidHeaderSpecific(X,VEC_CLASSID,3); 2146 if (!ts->ops->evaluatestep) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSEvaluateStep not implemented for type '%s'",((PetscObject)ts)->type_name); 2147 ierr = (*ts->ops->evaluatestep)(ts,order,X,done);CHKERRQ(ierr); 2148 PetscFunctionReturn(0); 2149 } 2150 2151 #undef __FUNCT__ 2152 #define __FUNCT__ "TSSolve" 2153 /*@ 2154 TSSolve - Steps the requested number of timesteps. 2155 2156 Collective on TS 2157 2158 Input Parameter: 2159 + ts - the TS context obtained from TSCreate() 2160 - x - the solution vector 2161 2162 Output Parameter: 2163 . ftime - time of the state vector x upon completion 2164 2165 Level: beginner 2166 2167 Notes: 2168 The final time returned by this function may be different from the time of the internally 2169 held state accessible by TSGetSolution() and TSGetTime() because the method may have 2170 stepped over the final time. 2171 2172 .keywords: TS, timestep, solve 2173 2174 .seealso: TSCreate(), TSSetSolution(), TSStep() 2175 @*/ 2176 PetscErrorCode TSSolve(TS ts,Vec x,PetscReal *ftime) 2177 { 2178 PetscBool flg; 2179 char filename[PETSC_MAX_PATH_LEN]; 2180 PetscViewer viewer; 2181 PetscErrorCode ierr; 2182 2183 PetscFunctionBegin; 2184 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2185 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 2186 if (ts->exact_final_time) { /* Need ts->vec_sol to be distinct so it is not overwritten when we interpolate at the end */ 2187 if (!ts->vec_sol || x == ts->vec_sol) { 2188 Vec y; 2189 ierr = VecDuplicate(x,&y);CHKERRQ(ierr); 2190 ierr = TSSetSolution(ts,y);CHKERRQ(ierr); 2191 ierr = VecDestroy(&y);CHKERRQ(ierr); /* grant ownership */ 2192 } 2193 ierr = VecCopy(x,ts->vec_sol);CHKERRQ(ierr); 2194 } else { 2195 ierr = TSSetSolution(ts,x);CHKERRQ(ierr); 2196 } 2197 ierr = TSSetUp(ts);CHKERRQ(ierr); 2198 /* reset time step and iteration counters */ 2199 ts->steps = 0; 2200 ts->ksp_its = 0; 2201 ts->snes_its = 0; 2202 ts->num_snes_failures = 0; 2203 ts->reject = 0; 2204 ts->reason = TS_CONVERGED_ITERATING; 2205 2206 if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */ 2207 ierr = (*ts->ops->solve)(ts);CHKERRQ(ierr); 2208 ierr = VecCopy(ts->vec_sol,x);CHKERRQ(ierr); 2209 if (ftime) *ftime = ts->ptime; 2210 } else { 2211 /* steps the requested number of timesteps. */ 2212 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 2213 if (ts->steps >= ts->max_steps) 2214 ts->reason = TS_CONVERGED_ITS; 2215 else if (ts->ptime >= ts->max_time) 2216 ts->reason = TS_CONVERGED_TIME; 2217 while (!ts->reason) { 2218 ierr = TSStep(ts);CHKERRQ(ierr); 2219 ierr = TSPostStep(ts);CHKERRQ(ierr); 2220 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 2221 } 2222 if (ts->exact_final_time && ts->ptime > ts->max_time) { 2223 ierr = TSInterpolate(ts,ts->max_time,x);CHKERRQ(ierr); 2224 if (ftime) *ftime = ts->max_time; 2225 } else { 2226 ierr = VecCopy(ts->vec_sol,x);CHKERRQ(ierr); 2227 if (ftime) *ftime = ts->ptime; 2228 } 2229 } 2230 ierr = PetscOptionsGetString(((PetscObject)ts)->prefix,"-ts_view",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 2231 if (flg && !PetscPreLoadingOn) { 2232 ierr = PetscViewerASCIIOpen(((PetscObject)ts)->comm,filename,&viewer);CHKERRQ(ierr); 2233 ierr = TSView(ts,viewer);CHKERRQ(ierr); 2234 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2235 } 2236 PetscFunctionReturn(0); 2237 } 2238 2239 #undef __FUNCT__ 2240 #define __FUNCT__ "TSMonitor" 2241 /*@ 2242 TSMonitor - Runs all user-provided monitor routines set using TSMonitorSet() 2243 2244 Collective on TS 2245 2246 Input Parameters: 2247 + ts - time stepping context obtained from TSCreate() 2248 . step - step number that has just completed 2249 . ptime - model time of the state 2250 - x - state at the current model time 2251 2252 Notes: 2253 TSMonitor() is typically used within the time stepping implementations. 2254 Users might call this function when using the TSStep() interface instead of TSSolve(). 2255 2256 Level: advanced 2257 2258 .keywords: TS, timestep 2259 @*/ 2260 PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec x) 2261 { 2262 PetscErrorCode ierr; 2263 PetscInt i,n = ts->numbermonitors; 2264 2265 PetscFunctionBegin; 2266 for (i=0; i<n; i++) { 2267 ierr = (*ts->monitor[i])(ts,step,ptime,x,ts->monitorcontext[i]);CHKERRQ(ierr); 2268 } 2269 PetscFunctionReturn(0); 2270 } 2271 2272 /* ------------------------------------------------------------------------*/ 2273 2274 #undef __FUNCT__ 2275 #define __FUNCT__ "TSMonitorLGCreate" 2276 /*@C 2277 TSMonitorLGCreate - Creates a line graph context for use with 2278 TS to monitor convergence of preconditioned residual norms. 2279 2280 Collective on TS 2281 2282 Input Parameters: 2283 + host - the X display to open, or null for the local machine 2284 . label - the title to put in the title bar 2285 . x, y - the screen coordinates of the upper left coordinate of the window 2286 - m, n - the screen width and height in pixels 2287 2288 Output Parameter: 2289 . draw - the drawing context 2290 2291 Options Database Key: 2292 . -ts_monitor_draw - automatically sets line graph monitor 2293 2294 Notes: 2295 Use TSMonitorLGDestroy() to destroy this line graph, not PetscDrawLGDestroy(). 2296 2297 Level: intermediate 2298 2299 .keywords: TS, monitor, line graph, residual, seealso 2300 2301 .seealso: TSMonitorLGDestroy(), TSMonitorSet() 2302 2303 @*/ 2304 PetscErrorCode TSMonitorLGCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 2305 { 2306 PetscDraw win; 2307 PetscErrorCode ierr; 2308 PetscDrawAxis axis; 2309 2310 PetscFunctionBegin; 2311 ierr = PetscDrawCreate(comm,host,label,x,y,m,n,&win);CHKERRQ(ierr); 2312 ierr = PetscDrawSetType(win,PETSC_DRAW_X);CHKERRQ(ierr); 2313 ierr = PetscDrawLGCreate(win,1,draw);CHKERRQ(ierr); 2314 ierr = PetscDrawLGIndicateDataPoints(*draw);CHKERRQ(ierr); 2315 ierr = PetscDrawLGGetAxis(*draw,&axis);CHKERRQ(ierr); 2316 ierr = PetscDrawAxisSetLabels(axis,"Integration Time","Iteration","Time");CHKERRQ(ierr); 2317 ierr = PetscLogObjectParent(*draw,win);CHKERRQ(ierr); 2318 PetscFunctionReturn(0); 2319 } 2320 2321 #undef __FUNCT__ 2322 #define __FUNCT__ "TSMonitorLG" 2323 PetscErrorCode TSMonitorLG(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 2324 { 2325 PetscDrawLG lg = (PetscDrawLG) monctx; 2326 PetscReal x,y = ptime; 2327 PetscErrorCode ierr; 2328 2329 PetscFunctionBegin; 2330 if (!monctx) { 2331 MPI_Comm comm; 2332 PetscViewer viewer; 2333 PetscDrawAxis axis; 2334 2335 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 2336 viewer = PETSC_VIEWER_DRAW_(comm); 2337 ierr = PetscViewerDrawGetDrawLG(viewer,0,&lg);CHKERRQ(ierr); 2338 ierr = PetscDrawLGGetAxis(lg,&axis);CHKERRQ(ierr); 2339 ierr = PetscDrawAxisSetLabels(axis,"Integration Time","Iteration","Time");CHKERRQ(ierr); 2340 } 2341 2342 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 2343 x = (PetscReal)n; 2344 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 2345 if (n < 20 || (n % 5)) { 2346 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 2347 } 2348 PetscFunctionReturn(0); 2349 } 2350 2351 #undef __FUNCT__ 2352 #define __FUNCT__ "TSMonitorLGDestroy" 2353 /*@C 2354 TSMonitorLGDestroy - Destroys a line graph context that was created 2355 with TSMonitorLGCreate(). 2356 2357 Collective on PetscDrawLG 2358 2359 Input Parameter: 2360 . draw - the drawing context 2361 2362 Level: intermediate 2363 2364 .keywords: TS, monitor, line graph, destroy 2365 2366 .seealso: TSMonitorLGCreate(), TSMonitorSet(), TSMonitorLG(); 2367 @*/ 2368 PetscErrorCode TSMonitorLGDestroy(PetscDrawLG *drawlg) 2369 { 2370 PetscDraw draw; 2371 PetscErrorCode ierr; 2372 2373 PetscFunctionBegin; 2374 ierr = PetscDrawLGGetDraw(*drawlg,&draw);CHKERRQ(ierr); 2375 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 2376 ierr = PetscDrawLGDestroy(drawlg);CHKERRQ(ierr); 2377 PetscFunctionReturn(0); 2378 } 2379 2380 #undef __FUNCT__ 2381 #define __FUNCT__ "TSGetTime" 2382 /*@ 2383 TSGetTime - Gets the time of the most recently completed step. 2384 2385 Not Collective 2386 2387 Input Parameter: 2388 . ts - the TS context obtained from TSCreate() 2389 2390 Output Parameter: 2391 . t - the current time 2392 2393 Level: beginner 2394 2395 Note: 2396 When called during time step evaluation (e.g. during residual evaluation or via hooks set using TSSetPreStep(), 2397 TSSetPreStage(), or TSSetPostStep()), the time is the time at the start of the step being evaluated. 2398 2399 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 2400 2401 .keywords: TS, get, time 2402 @*/ 2403 PetscErrorCode TSGetTime(TS ts,PetscReal* t) 2404 { 2405 PetscFunctionBegin; 2406 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2407 PetscValidRealPointer(t,2); 2408 *t = ts->ptime; 2409 PetscFunctionReturn(0); 2410 } 2411 2412 #undef __FUNCT__ 2413 #define __FUNCT__ "TSSetTime" 2414 /*@ 2415 TSSetTime - Allows one to reset the time. 2416 2417 Logically Collective on TS 2418 2419 Input Parameters: 2420 + ts - the TS context obtained from TSCreate() 2421 - time - the time 2422 2423 Level: intermediate 2424 2425 .seealso: TSGetTime(), TSSetDuration() 2426 2427 .keywords: TS, set, time 2428 @*/ 2429 PetscErrorCode TSSetTime(TS ts, PetscReal t) 2430 { 2431 PetscFunctionBegin; 2432 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2433 PetscValidLogicalCollectiveReal(ts,t,2); 2434 ts->ptime = t; 2435 PetscFunctionReturn(0); 2436 } 2437 2438 #undef __FUNCT__ 2439 #define __FUNCT__ "TSSetOptionsPrefix" 2440 /*@C 2441 TSSetOptionsPrefix - Sets the prefix used for searching for all 2442 TS options in the database. 2443 2444 Logically Collective on TS 2445 2446 Input Parameter: 2447 + ts - The TS context 2448 - prefix - The prefix to prepend to all option names 2449 2450 Notes: 2451 A hyphen (-) must NOT be given at the beginning of the prefix name. 2452 The first character of all runtime options is AUTOMATICALLY the 2453 hyphen. 2454 2455 Level: advanced 2456 2457 .keywords: TS, set, options, prefix, database 2458 2459 .seealso: TSSetFromOptions() 2460 2461 @*/ 2462 PetscErrorCode TSSetOptionsPrefix(TS ts,const char prefix[]) 2463 { 2464 PetscErrorCode ierr; 2465 SNES snes; 2466 2467 PetscFunctionBegin; 2468 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2469 ierr = PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2470 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2471 ierr = SNESSetOptionsPrefix(snes,prefix);CHKERRQ(ierr); 2472 PetscFunctionReturn(0); 2473 } 2474 2475 2476 #undef __FUNCT__ 2477 #define __FUNCT__ "TSAppendOptionsPrefix" 2478 /*@C 2479 TSAppendOptionsPrefix - Appends to the prefix used for searching for all 2480 TS options in the database. 2481 2482 Logically Collective on TS 2483 2484 Input Parameter: 2485 + ts - The TS context 2486 - prefix - The prefix to prepend to all option names 2487 2488 Notes: 2489 A hyphen (-) must NOT be given at the beginning of the prefix name. 2490 The first character of all runtime options is AUTOMATICALLY the 2491 hyphen. 2492 2493 Level: advanced 2494 2495 .keywords: TS, append, options, prefix, database 2496 2497 .seealso: TSGetOptionsPrefix() 2498 2499 @*/ 2500 PetscErrorCode TSAppendOptionsPrefix(TS ts,const char prefix[]) 2501 { 2502 PetscErrorCode ierr; 2503 SNES snes; 2504 2505 PetscFunctionBegin; 2506 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2507 ierr = PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2508 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2509 ierr = SNESAppendOptionsPrefix(snes,prefix);CHKERRQ(ierr); 2510 PetscFunctionReturn(0); 2511 } 2512 2513 #undef __FUNCT__ 2514 #define __FUNCT__ "TSGetOptionsPrefix" 2515 /*@C 2516 TSGetOptionsPrefix - Sets the prefix used for searching for all 2517 TS options in the database. 2518 2519 Not Collective 2520 2521 Input Parameter: 2522 . ts - The TS context 2523 2524 Output Parameter: 2525 . prefix - A pointer to the prefix string used 2526 2527 Notes: On the fortran side, the user should pass in a string 'prifix' of 2528 sufficient length to hold the prefix. 2529 2530 Level: intermediate 2531 2532 .keywords: TS, get, options, prefix, database 2533 2534 .seealso: TSAppendOptionsPrefix() 2535 @*/ 2536 PetscErrorCode TSGetOptionsPrefix(TS ts,const char *prefix[]) 2537 { 2538 PetscErrorCode ierr; 2539 2540 PetscFunctionBegin; 2541 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2542 PetscValidPointer(prefix,2); 2543 ierr = PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2544 PetscFunctionReturn(0); 2545 } 2546 2547 #undef __FUNCT__ 2548 #define __FUNCT__ "TSGetRHSJacobian" 2549 /*@C 2550 TSGetRHSJacobian - Returns the Jacobian J at the present timestep. 2551 2552 Not Collective, but parallel objects are returned if TS is parallel 2553 2554 Input Parameter: 2555 . ts - The TS context obtained from TSCreate() 2556 2557 Output Parameters: 2558 + J - The Jacobian J of F, where U_t = F(U,t) 2559 . M - The preconditioner matrix, usually the same as J 2560 . func - Function to compute the Jacobian of the RHS 2561 - ctx - User-defined context for Jacobian evaluation routine 2562 2563 Notes: You can pass in PETSC_NULL for any return argument you do not need. 2564 2565 Level: intermediate 2566 2567 .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 2568 2569 .keywords: TS, timestep, get, matrix, Jacobian 2570 @*/ 2571 PetscErrorCode TSGetRHSJacobian(TS ts,Mat *J,Mat *M,TSRHSJacobian *func,void **ctx) 2572 { 2573 PetscErrorCode ierr; 2574 SNES snes; 2575 DM dm; 2576 2577 PetscFunctionBegin; 2578 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2579 ierr = SNESGetJacobian(snes,J,M,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2580 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 2581 ierr = DMTSGetRHSJacobian(dm,func,ctx);CHKERRQ(ierr); 2582 PetscFunctionReturn(0); 2583 } 2584 2585 #undef __FUNCT__ 2586 #define __FUNCT__ "TSGetIJacobian" 2587 /*@C 2588 TSGetIJacobian - Returns the implicit Jacobian at the present timestep. 2589 2590 Not Collective, but parallel objects are returned if TS is parallel 2591 2592 Input Parameter: 2593 . ts - The TS context obtained from TSCreate() 2594 2595 Output Parameters: 2596 + A - The Jacobian of F(t,U,U_t) 2597 . B - The preconditioner matrix, often the same as A 2598 . f - The function to compute the matrices 2599 - ctx - User-defined context for Jacobian evaluation routine 2600 2601 Notes: You can pass in PETSC_NULL for any return argument you do not need. 2602 2603 Level: advanced 2604 2605 .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 2606 2607 .keywords: TS, timestep, get, matrix, Jacobian 2608 @*/ 2609 PetscErrorCode TSGetIJacobian(TS ts,Mat *A,Mat *B,TSIJacobian *f,void **ctx) 2610 { 2611 PetscErrorCode ierr; 2612 SNES snes; 2613 DM dm; 2614 PetscFunctionBegin; 2615 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2616 ierr = SNESGetJacobian(snes,A,B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2617 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 2618 ierr = DMTSGetIJacobian(dm,f,ctx);CHKERRQ(ierr); 2619 PetscFunctionReturn(0); 2620 } 2621 2622 typedef struct { 2623 PetscViewer viewer; 2624 Vec initialsolution; 2625 PetscBool showinitial; 2626 } TSMonitorSolutionCtx; 2627 2628 #undef __FUNCT__ 2629 #define __FUNCT__ "TSMonitorSolution" 2630 /*@C 2631 TSMonitorSolution - Monitors progress of the TS solvers by calling 2632 VecView() for the solution at each timestep 2633 2634 Collective on TS 2635 2636 Input Parameters: 2637 + ts - the TS context 2638 . step - current time-step 2639 . ptime - current time 2640 - dummy - either a viewer or PETSC_NULL 2641 2642 Level: intermediate 2643 2644 .keywords: TS, vector, monitor, view 2645 2646 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 2647 @*/ 2648 PetscErrorCode TSMonitorSolution(TS ts,PetscInt step,PetscReal ptime,Vec x,void *dummy) 2649 { 2650 PetscErrorCode ierr; 2651 TSMonitorSolutionCtx *ictx = (TSMonitorSolutionCtx*)dummy; 2652 2653 PetscFunctionBegin; 2654 if (!step && ictx->showinitial) { 2655 if (!ictx->initialsolution) { 2656 ierr = VecDuplicate(x,&ictx->initialsolution);CHKERRQ(ierr); 2657 } 2658 ierr = VecCopy(x,ictx->initialsolution);CHKERRQ(ierr); 2659 } 2660 if (ictx->showinitial) { 2661 PetscReal pause; 2662 ierr = PetscViewerDrawGetPause(ictx->viewer,&pause);CHKERRQ(ierr); 2663 ierr = PetscViewerDrawSetPause(ictx->viewer,0.0);CHKERRQ(ierr); 2664 ierr = VecView(ictx->initialsolution,ictx->viewer);CHKERRQ(ierr); 2665 ierr = PetscViewerDrawSetPause(ictx->viewer,pause);CHKERRQ(ierr); 2666 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);CHKERRQ(ierr); 2667 } 2668 ierr = VecView(x,ictx->viewer);CHKERRQ(ierr); 2669 if (ictx->showinitial) { 2670 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);CHKERRQ(ierr); 2671 } 2672 PetscFunctionReturn(0); 2673 } 2674 2675 2676 #undef __FUNCT__ 2677 #define __FUNCT__ "TSMonitorSolutionDestroy" 2678 /*@C 2679 TSMonitorSolutionDestroy - Destroys the monitor context for TSMonitorSolution 2680 2681 Collective on TS 2682 2683 Input Parameters: 2684 . ctx - the monitor context 2685 2686 Level: intermediate 2687 2688 .keywords: TS, vector, monitor, view 2689 2690 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorSolution() 2691 @*/ 2692 PetscErrorCode TSMonitorSolutionDestroy(void **ctx) 2693 { 2694 PetscErrorCode ierr; 2695 TSMonitorSolutionCtx *ictx = *(TSMonitorSolutionCtx**)ctx; 2696 2697 PetscFunctionBegin; 2698 ierr = PetscViewerDestroy(&ictx->viewer);CHKERRQ(ierr); 2699 ierr = VecDestroy(&ictx->initialsolution);CHKERRQ(ierr); 2700 ierr = PetscFree(ictx);CHKERRQ(ierr); 2701 PetscFunctionReturn(0); 2702 } 2703 2704 #undef __FUNCT__ 2705 #define __FUNCT__ "TSMonitorSolutionCreate" 2706 /*@C 2707 TSMonitorSolutionCreate - Creates the monitor context for TSMonitorSolution 2708 2709 Collective on TS 2710 2711 Input Parameter: 2712 . ts - time-step context 2713 2714 Output Patameter: 2715 . ctx - the monitor context 2716 2717 Level: intermediate 2718 2719 .keywords: TS, vector, monitor, view 2720 2721 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorSolution() 2722 @*/ 2723 PetscErrorCode TSMonitorSolutionCreate(TS ts,PetscViewer viewer,void **ctx) 2724 { 2725 PetscErrorCode ierr; 2726 TSMonitorSolutionCtx *ictx; 2727 2728 PetscFunctionBegin; 2729 ierr = PetscNew(TSMonitorSolutionCtx,&ictx);CHKERRQ(ierr); 2730 *ctx = (void*)ictx; 2731 if (!viewer) { 2732 viewer = PETSC_VIEWER_DRAW_(((PetscObject)ts)->comm); 2733 } 2734 ierr = PetscObjectReference((PetscObject)viewer);CHKERRQ(ierr); 2735 ictx->viewer = viewer; 2736 ictx->showinitial = PETSC_FALSE; 2737 ierr = PetscOptionsGetBool(((PetscObject)ts)->prefix,"-ts_monitor_solution_initial",&ictx->showinitial,PETSC_NULL);CHKERRQ(ierr); 2738 PetscFunctionReturn(0); 2739 } 2740 2741 #undef __FUNCT__ 2742 #define __FUNCT__ "TSSetDM" 2743 /*@ 2744 TSSetDM - Sets the DM that may be used by some preconditioners 2745 2746 Logically Collective on TS and DM 2747 2748 Input Parameters: 2749 + ts - the preconditioner context 2750 - dm - the dm 2751 2752 Level: intermediate 2753 2754 2755 .seealso: TSGetDM(), SNESSetDM(), SNESGetDM() 2756 @*/ 2757 PetscErrorCode TSSetDM(TS ts,DM dm) 2758 { 2759 PetscErrorCode ierr; 2760 SNES snes; 2761 TSDM tsdm; 2762 2763 PetscFunctionBegin; 2764 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2765 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 2766 if (ts->dm) { /* Move the TSDM context over to the new DM unless the new DM already has one */ 2767 PetscContainer oldcontainer,container; 2768 ierr = PetscObjectQuery((PetscObject)ts->dm,"TSDM",(PetscObject*)&oldcontainer);CHKERRQ(ierr); 2769 ierr = PetscObjectQuery((PetscObject)dm,"TSDM",(PetscObject*)&container);CHKERRQ(ierr); 2770 if (oldcontainer && !container) { 2771 ierr = DMTSCopyContext(ts->dm,dm);CHKERRQ(ierr); 2772 ierr = DMTSGetContext(ts->dm,&tsdm);CHKERRQ(ierr); 2773 if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */ 2774 tsdm->originaldm = dm; 2775 } 2776 } 2777 ierr = DMDestroy(&ts->dm);CHKERRQ(ierr); 2778 } 2779 ts->dm = dm; 2780 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2781 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 2782 PetscFunctionReturn(0); 2783 } 2784 2785 #undef __FUNCT__ 2786 #define __FUNCT__ "TSGetDM" 2787 /*@ 2788 TSGetDM - Gets the DM that may be used by some preconditioners 2789 2790 Not Collective 2791 2792 Input Parameter: 2793 . ts - the preconditioner context 2794 2795 Output Parameter: 2796 . dm - the dm 2797 2798 Level: intermediate 2799 2800 2801 .seealso: TSSetDM(), SNESSetDM(), SNESGetDM() 2802 @*/ 2803 PetscErrorCode TSGetDM(TS ts,DM *dm) 2804 { 2805 PetscErrorCode ierr; 2806 2807 PetscFunctionBegin; 2808 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2809 if (!ts->dm) { 2810 ierr = DMShellCreate(((PetscObject)ts)->comm,&ts->dm);CHKERRQ(ierr); 2811 if (ts->snes) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 2812 } 2813 *dm = ts->dm; 2814 PetscFunctionReturn(0); 2815 } 2816 2817 #undef __FUNCT__ 2818 #define __FUNCT__ "SNESTSFormFunction" 2819 /*@ 2820 SNESTSFormFunction - Function to evaluate nonlinear residual 2821 2822 Logically Collective on SNES 2823 2824 Input Parameter: 2825 + snes - nonlinear solver 2826 . X - the current state at which to evaluate the residual 2827 - ctx - user context, must be a TS 2828 2829 Output Parameter: 2830 . F - the nonlinear residual 2831 2832 Notes: 2833 This function is not normally called by users and is automatically registered with the SNES used by TS. 2834 It is most frequently passed to MatFDColoringSetFunction(). 2835 2836 Level: advanced 2837 2838 .seealso: SNESSetFunction(), MatFDColoringSetFunction() 2839 @*/ 2840 PetscErrorCode SNESTSFormFunction(SNES snes,Vec X,Vec F,void *ctx) 2841 { 2842 TS ts = (TS)ctx; 2843 PetscErrorCode ierr; 2844 2845 PetscFunctionBegin; 2846 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2847 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2848 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 2849 PetscValidHeaderSpecific(ts,TS_CLASSID,4); 2850 ierr = (ts->ops->snesfunction)(snes,X,F,ts);CHKERRQ(ierr); 2851 PetscFunctionReturn(0); 2852 } 2853 2854 #undef __FUNCT__ 2855 #define __FUNCT__ "SNESTSFormJacobian" 2856 /*@ 2857 SNESTSFormJacobian - Function to evaluate the Jacobian 2858 2859 Collective on SNES 2860 2861 Input Parameter: 2862 + snes - nonlinear solver 2863 . X - the current state at which to evaluate the residual 2864 - ctx - user context, must be a TS 2865 2866 Output Parameter: 2867 + A - the Jacobian 2868 . B - the preconditioning matrix (may be the same as A) 2869 - flag - indicates any structure change in the matrix 2870 2871 Notes: 2872 This function is not normally called by users and is automatically registered with the SNES used by TS. 2873 2874 Level: developer 2875 2876 .seealso: SNESSetJacobian() 2877 @*/ 2878 PetscErrorCode SNESTSFormJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flag,void *ctx) 2879 { 2880 TS ts = (TS)ctx; 2881 PetscErrorCode ierr; 2882 2883 PetscFunctionBegin; 2884 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2885 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2886 PetscValidPointer(A,3); 2887 PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 2888 PetscValidPointer(B,4); 2889 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); 2890 PetscValidPointer(flag,5); 2891 PetscValidHeaderSpecific(ts,TS_CLASSID,6); 2892 ierr = (ts->ops->snesjacobian)(snes,X,A,B,flag,ts);CHKERRQ(ierr); 2893 PetscFunctionReturn(0); 2894 } 2895 2896 #undef __FUNCT__ 2897 #define __FUNCT__ "TSComputeRHSFunctionLinear" 2898 /*@C 2899 TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems only 2900 2901 Collective on TS 2902 2903 Input Arguments: 2904 + ts - time stepping context 2905 . t - time at which to evaluate 2906 . X - state at which to evaluate 2907 - ctx - context 2908 2909 Output Arguments: 2910 . F - right hand side 2911 2912 Level: intermediate 2913 2914 Notes: 2915 This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems. 2916 The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian(). 2917 2918 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant() 2919 @*/ 2920 PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec X,Vec F,void *ctx) 2921 { 2922 PetscErrorCode ierr; 2923 Mat Arhs,Brhs; 2924 MatStructure flg2; 2925 2926 PetscFunctionBegin; 2927 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 2928 ierr = TSComputeRHSJacobian(ts,t,X,&Arhs,&Brhs,&flg2);CHKERRQ(ierr); 2929 ierr = MatMult(Arhs,X,F);CHKERRQ(ierr); 2930 PetscFunctionReturn(0); 2931 } 2932 2933 #undef __FUNCT__ 2934 #define __FUNCT__ "TSComputeRHSJacobianConstant" 2935 /*@C 2936 TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent. 2937 2938 Collective on TS 2939 2940 Input Arguments: 2941 + ts - time stepping context 2942 . t - time at which to evaluate 2943 . X - state at which to evaluate 2944 - ctx - context 2945 2946 Output Arguments: 2947 + A - pointer to operator 2948 . B - pointer to preconditioning matrix 2949 - flg - matrix structure flag 2950 2951 Level: intermediate 2952 2953 Notes: 2954 This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems. 2955 2956 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear() 2957 @*/ 2958 PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec X,Mat *A,Mat *B,MatStructure *flg,void *ctx) 2959 { 2960 2961 PetscFunctionBegin; 2962 *flg = SAME_PRECONDITIONER; 2963 PetscFunctionReturn(0); 2964 } 2965 2966 #undef __FUNCT__ 2967 #define __FUNCT__ "TSComputeIFunctionLinear" 2968 /*@C 2969 TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only 2970 2971 Collective on TS 2972 2973 Input Arguments: 2974 + ts - time stepping context 2975 . t - time at which to evaluate 2976 . X - state at which to evaluate 2977 . Xdot - time derivative of state vector 2978 - ctx - context 2979 2980 Output Arguments: 2981 . F - left hand side 2982 2983 Level: intermediate 2984 2985 Notes: 2986 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 2987 user is required to write their own TSComputeIFunction. 2988 This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems. 2989 The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian(). 2990 2991 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant() 2992 @*/ 2993 PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx) 2994 { 2995 PetscErrorCode ierr; 2996 Mat A,B; 2997 MatStructure flg2; 2998 2999 PetscFunctionBegin; 3000 ierr = TSGetIJacobian(ts,&A,&B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 3001 ierr = TSComputeIJacobian(ts,t,X,Xdot,1.0,&A,&B,&flg2,PETSC_TRUE);CHKERRQ(ierr); 3002 ierr = MatMult(A,Xdot,F);CHKERRQ(ierr); 3003 PetscFunctionReturn(0); 3004 } 3005 3006 #undef __FUNCT__ 3007 #define __FUNCT__ "TSComputeIJacobianConstant" 3008 /*@C 3009 TSComputeIJacobianConstant - Reuses a Jacobian that is time-independent. 3010 3011 Collective on TS 3012 3013 Input Arguments: 3014 + ts - time stepping context 3015 . t - time at which to evaluate 3016 . X - state at which to evaluate 3017 . Xdot - time derivative of state vector 3018 . shift - shift to apply 3019 - ctx - context 3020 3021 Output Arguments: 3022 + A - pointer to operator 3023 . B - pointer to preconditioning matrix 3024 - flg - matrix structure flag 3025 3026 Level: intermediate 3027 3028 Notes: 3029 This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems. 3030 3031 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear() 3032 @*/ 3033 PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,void *ctx) 3034 { 3035 3036 PetscFunctionBegin; 3037 *flg = SAME_PRECONDITIONER; 3038 PetscFunctionReturn(0); 3039 } 3040 3041 3042 #undef __FUNCT__ 3043 #define __FUNCT__ "TSGetConvergedReason" 3044 /*@ 3045 TSGetConvergedReason - Gets the reason the TS iteration was stopped. 3046 3047 Not Collective 3048 3049 Input Parameter: 3050 . ts - the TS context 3051 3052 Output Parameter: 3053 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 3054 manual pages for the individual convergence tests for complete lists 3055 3056 Level: intermediate 3057 3058 Notes: 3059 Can only be called after the call to TSSolve() is complete. 3060 3061 .keywords: TS, nonlinear, set, convergence, test 3062 3063 .seealso: TSSetConvergenceTest(), TSConvergedReason 3064 @*/ 3065 PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason) 3066 { 3067 PetscFunctionBegin; 3068 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3069 PetscValidPointer(reason,2); 3070 *reason = ts->reason; 3071 PetscFunctionReturn(0); 3072 } 3073 3074 #undef __FUNCT__ 3075 #define __FUNCT__ "TSGetSNESIterations" 3076 /*@ 3077 TSGetSNESIterations - Gets the total number of nonlinear iterations 3078 used by the time integrator. 3079 3080 Not Collective 3081 3082 Input Parameter: 3083 . ts - TS context 3084 3085 Output Parameter: 3086 . nits - number of nonlinear iterations 3087 3088 Notes: 3089 This counter is reset to zero for each successive call to TSSolve(). 3090 3091 Level: intermediate 3092 3093 .keywords: TS, get, number, nonlinear, iterations 3094 3095 .seealso: TSGetKSPIterations() 3096 @*/ 3097 PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits) 3098 { 3099 PetscFunctionBegin; 3100 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3101 PetscValidIntPointer(nits,2); 3102 *nits = ts->snes_its; 3103 PetscFunctionReturn(0); 3104 } 3105 3106 #undef __FUNCT__ 3107 #define __FUNCT__ "TSGetKSPIterations" 3108 /*@ 3109 TSGetKSPIterations - Gets the total number of linear iterations 3110 used by the time integrator. 3111 3112 Not Collective 3113 3114 Input Parameter: 3115 . ts - TS context 3116 3117 Output Parameter: 3118 . lits - number of linear iterations 3119 3120 Notes: 3121 This counter is reset to zero for each successive call to TSSolve(). 3122 3123 Level: intermediate 3124 3125 .keywords: TS, get, number, linear, iterations 3126 3127 .seealso: TSGetSNESIterations(), SNESGetKSPIterations() 3128 @*/ 3129 PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits) 3130 { 3131 PetscFunctionBegin; 3132 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3133 PetscValidIntPointer(lits,2); 3134 *lits = ts->ksp_its; 3135 PetscFunctionReturn(0); 3136 } 3137 3138 #undef __FUNCT__ 3139 #define __FUNCT__ "TSGetStepRejections" 3140 /*@ 3141 TSGetStepRejections - Gets the total number of rejected steps. 3142 3143 Not Collective 3144 3145 Input Parameter: 3146 . ts - TS context 3147 3148 Output Parameter: 3149 . rejects - number of steps rejected 3150 3151 Notes: 3152 This counter is reset to zero for each successive call to TSSolve(). 3153 3154 Level: intermediate 3155 3156 .keywords: TS, get, number 3157 3158 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails() 3159 @*/ 3160 PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects) 3161 { 3162 PetscFunctionBegin; 3163 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3164 PetscValidIntPointer(rejects,2); 3165 *rejects = ts->reject; 3166 PetscFunctionReturn(0); 3167 } 3168 3169 #undef __FUNCT__ 3170 #define __FUNCT__ "TSGetSNESFailures" 3171 /*@ 3172 TSGetSNESFailures - Gets the total number of failed SNES solves 3173 3174 Not Collective 3175 3176 Input Parameter: 3177 . ts - TS context 3178 3179 Output Parameter: 3180 . fails - number of failed nonlinear solves 3181 3182 Notes: 3183 This counter is reset to zero for each successive call to TSSolve(). 3184 3185 Level: intermediate 3186 3187 .keywords: TS, get, number 3188 3189 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures() 3190 @*/ 3191 PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails) 3192 { 3193 PetscFunctionBegin; 3194 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3195 PetscValidIntPointer(fails,2); 3196 *fails = ts->num_snes_failures; 3197 PetscFunctionReturn(0); 3198 } 3199 3200 #undef __FUNCT__ 3201 #define __FUNCT__ "TSSetMaxStepRejections" 3202 /*@ 3203 TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails 3204 3205 Not Collective 3206 3207 Input Parameter: 3208 + ts - TS context 3209 - rejects - maximum number of rejected steps, pass -1 for unlimited 3210 3211 Notes: 3212 The counter is reset to zero for each step 3213 3214 Options Database Key: 3215 . -ts_max_reject - Maximum number of step rejections before a step fails 3216 3217 Level: intermediate 3218 3219 .keywords: TS, set, maximum, number 3220 3221 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 3222 @*/ 3223 PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects) 3224 { 3225 PetscFunctionBegin; 3226 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3227 ts->max_reject = rejects; 3228 PetscFunctionReturn(0); 3229 } 3230 3231 #undef __FUNCT__ 3232 #define __FUNCT__ "TSSetMaxSNESFailures" 3233 /*@ 3234 TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves 3235 3236 Not Collective 3237 3238 Input Parameter: 3239 + ts - TS context 3240 - fails - maximum number of failed nonlinear solves, pass -1 for unlimited 3241 3242 Notes: 3243 The counter is reset to zero for each successive call to TSSolve(). 3244 3245 Options Database Key: 3246 . -ts_max_snes_failures - Maximum number of nonlinear solve failures 3247 3248 Level: intermediate 3249 3250 .keywords: TS, set, maximum, number 3251 3252 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason() 3253 @*/ 3254 PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails) 3255 { 3256 PetscFunctionBegin; 3257 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3258 ts->max_snes_failures = fails; 3259 PetscFunctionReturn(0); 3260 } 3261 3262 #undef __FUNCT__ 3263 #define __FUNCT__ "TSSetErrorIfStepFails()" 3264 /*@ 3265 TSSetErrorIfStepFails - Error if no step succeeds 3266 3267 Not Collective 3268 3269 Input Parameter: 3270 + ts - TS context 3271 - err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure 3272 3273 Options Database Key: 3274 . -ts_error_if_step_fails - Error if no step succeeds 3275 3276 Level: intermediate 3277 3278 .keywords: TS, set, error 3279 3280 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 3281 @*/ 3282 PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err) 3283 { 3284 PetscFunctionBegin; 3285 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3286 ts->errorifstepfailed = err; 3287 PetscFunctionReturn(0); 3288 } 3289 3290 #undef __FUNCT__ 3291 #define __FUNCT__ "TSMonitorSolutionBinary" 3292 /*@C 3293 TSMonitorSolutionBinary - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file 3294 3295 Collective on TS 3296 3297 Input Parameters: 3298 + ts - the TS context 3299 . step - current time-step 3300 . ptime - current time 3301 . x - current state 3302 - viewer - binary viewer 3303 3304 Level: intermediate 3305 3306 .keywords: TS, vector, monitor, view 3307 3308 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3309 @*/ 3310 PetscErrorCode TSMonitorSolutionBinary(TS ts,PetscInt step,PetscReal ptime,Vec x,void *viewer) 3311 { 3312 PetscErrorCode ierr; 3313 PetscViewer v = (PetscViewer)viewer; 3314 3315 PetscFunctionBegin; 3316 ierr = VecView(x,v);CHKERRQ(ierr); 3317 PetscFunctionReturn(0); 3318 } 3319 3320 #undef __FUNCT__ 3321 #define __FUNCT__ "TSMonitorSolutionVTK" 3322 /*@C 3323 TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep. 3324 3325 Collective on TS 3326 3327 Input Parameters: 3328 + ts - the TS context 3329 . step - current time-step 3330 . ptime - current time 3331 . x - current state 3332 - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 3333 3334 Level: intermediate 3335 3336 Notes: 3337 The VTK format does not allow writing multiple time steps in the same file, therefore a different file will be written for each time step. 3338 These are named according to the file name template. 3339 3340 This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy(). 3341 3342 .keywords: TS, vector, monitor, view 3343 3344 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3345 @*/ 3346 PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec x,void *filenametemplate) 3347 { 3348 PetscErrorCode ierr; 3349 char filename[PETSC_MAX_PATH_LEN]; 3350 PetscViewer viewer; 3351 3352 PetscFunctionBegin; 3353 ierr = PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);CHKERRQ(ierr); 3354 ierr = PetscViewerVTKOpen(((PetscObject)ts)->comm,filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 3355 ierr = VecView(x,viewer);CHKERRQ(ierr); 3356 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3357 PetscFunctionReturn(0); 3358 } 3359 3360 #undef __FUNCT__ 3361 #define __FUNCT__ "TSMonitorSolutionVTKDestroy" 3362 /*@C 3363 TSMonitorSolutionVTKDestroy - Destroy context for monitoring 3364 3365 Collective on TS 3366 3367 Input Parameters: 3368 . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 3369 3370 Level: intermediate 3371 3372 Note: 3373 This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK(). 3374 3375 .keywords: TS, vector, monitor, view 3376 3377 .seealso: TSMonitorSet(), TSMonitorSolutionVTK() 3378 @*/ 3379 PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate) 3380 { 3381 PetscErrorCode ierr; 3382 3383 PetscFunctionBegin; 3384 ierr = PetscFree(*(char**)filenametemplate);CHKERRQ(ierr); 3385 PetscFunctionReturn(0); 3386 } 3387 3388 #undef __FUNCT__ 3389 #define __FUNCT__ "TSGetAdapt" 3390 /*@ 3391 TSGetAdapt - Get the adaptive controller context for the current method 3392 3393 Collective on TS if controller has not been created yet 3394 3395 Input Arguments: 3396 . ts - time stepping context 3397 3398 Output Arguments: 3399 . adapt - adaptive controller 3400 3401 Level: intermediate 3402 3403 .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose() 3404 @*/ 3405 PetscErrorCode TSGetAdapt(TS ts,TSAdapt *adapt) 3406 { 3407 PetscErrorCode ierr; 3408 3409 PetscFunctionBegin; 3410 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3411 PetscValidPointer(adapt,2); 3412 if (!ts->adapt) { 3413 ierr = TSAdaptCreate(((PetscObject)ts)->comm,&ts->adapt);CHKERRQ(ierr); 3414 ierr = PetscLogObjectParent(ts,ts->adapt);CHKERRQ(ierr); 3415 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);CHKERRQ(ierr); 3416 } 3417 *adapt = ts->adapt; 3418 PetscFunctionReturn(0); 3419 } 3420 3421 #undef __FUNCT__ 3422 #define __FUNCT__ "TSSetTolerances" 3423 /*@ 3424 TSSetTolerances - Set tolerances for local truncation error when using adaptive controller 3425 3426 Logically Collective 3427 3428 Input Arguments: 3429 + ts - time integration context 3430 . atol - scalar absolute tolerances, PETSC_DECIDE to leave current value 3431 . vatol - vector of absolute tolerances or PETSC_NULL, used in preference to atol if present 3432 . rtol - scalar relative tolerances, PETSC_DECIDE to leave current value 3433 - vrtol - vector of relative tolerances or PETSC_NULL, used in preference to atol if present 3434 3435 Level: beginner 3436 3437 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSGetTolerances() 3438 @*/ 3439 PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol) 3440 { 3441 PetscErrorCode ierr; 3442 3443 PetscFunctionBegin; 3444 if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol; 3445 if (vatol) { 3446 ierr = PetscObjectReference((PetscObject)vatol);CHKERRQ(ierr); 3447 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 3448 ts->vatol = vatol; 3449 } 3450 if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol; 3451 if (vrtol) { 3452 ierr = PetscObjectReference((PetscObject)vrtol);CHKERRQ(ierr); 3453 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 3454 ts->vrtol = vrtol; 3455 } 3456 PetscFunctionReturn(0); 3457 } 3458 3459 #undef __FUNCT__ 3460 #define __FUNCT__ "TSGetTolerances" 3461 /*@ 3462 TSGetTolerances - Get tolerances for local truncation error when using adaptive controller 3463 3464 Logically Collective 3465 3466 Input Arguments: 3467 . ts - time integration context 3468 3469 Output Arguments: 3470 + atol - scalar absolute tolerances, PETSC_NULL to ignore 3471 . vatol - vector of absolute tolerances, PETSC_NULL to ignore 3472 . rtol - scalar relative tolerances, PETSC_NULL to ignore 3473 - vrtol - vector of relative tolerances, PETSC_NULL to ignore 3474 3475 Level: beginner 3476 3477 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSSetTolerances() 3478 @*/ 3479 PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol) 3480 { 3481 3482 PetscFunctionBegin; 3483 if (atol) *atol = ts->atol; 3484 if (vatol) *vatol = ts->vatol; 3485 if (rtol) *rtol = ts->rtol; 3486 if (vrtol) *vrtol = ts->vrtol; 3487 PetscFunctionReturn(0); 3488 } 3489 3490 #undef __FUNCT__ 3491 #define __FUNCT__ "TSErrorNormWRMS" 3492 /*@ 3493 TSErrorNormWRMS - compute a weighted norm of the difference between a vector and the current state 3494 3495 Collective on TS 3496 3497 Input Arguments: 3498 + ts - time stepping context 3499 - Y - state vector to be compared to ts->vec_sol 3500 3501 Output Arguments: 3502 . norm - weighted norm, a value of 1.0 is considered small 3503 3504 Level: developer 3505 3506 .seealso: TSSetTolerances() 3507 @*/ 3508 PetscErrorCode TSErrorNormWRMS(TS ts,Vec Y,PetscReal *norm) 3509 { 3510 PetscErrorCode ierr; 3511 PetscInt i,n,N; 3512 const PetscScalar *x,*y; 3513 Vec X; 3514 PetscReal sum,gsum; 3515 3516 PetscFunctionBegin; 3517 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3518 PetscValidHeaderSpecific(Y,VEC_CLASSID,2); 3519 PetscValidPointer(norm,3); 3520 X = ts->vec_sol; 3521 PetscCheckSameTypeAndComm(X,1,Y,2); 3522 if (X == Y) SETERRQ(((PetscObject)X)->comm,PETSC_ERR_ARG_IDN,"Y cannot be the TS solution vector"); 3523 3524 ierr = VecGetSize(X,&N);CHKERRQ(ierr); 3525 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr); 3526 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 3527 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 3528 sum = 0.; 3529 if (ts->vatol && ts->vrtol) { 3530 const PetscScalar *atol,*rtol; 3531 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3532 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3533 for (i=0; i<n; i++) { 3534 PetscReal tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(x[i]),PetscAbsScalar(y[i])); 3535 sum += PetscSqr(PetscAbsScalar(y[i] - x[i]) / tol); 3536 } 3537 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3538 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3539 } else if (ts->vatol) { /* vector atol, scalar rtol */ 3540 const PetscScalar *atol; 3541 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3542 for (i=0; i<n; i++) { 3543 PetscReal tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(x[i]),PetscAbsScalar(y[i])); 3544 sum += PetscSqr(PetscAbsScalar(y[i] - x[i]) / tol); 3545 } 3546 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3547 } else if (ts->vrtol) { /* scalar atol, vector rtol */ 3548 const PetscScalar *rtol; 3549 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3550 for (i=0; i<n; i++) { 3551 PetscReal tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(x[i]),PetscAbsScalar(y[i])); 3552 sum += PetscSqr(PetscAbsScalar(y[i] - x[i]) / tol); 3553 } 3554 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3555 } else { /* scalar atol, scalar rtol */ 3556 for (i=0; i<n; i++) { 3557 PetscReal tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(x[i]),PetscAbsScalar(y[i])); 3558 sum += PetscSqr(PetscAbsScalar(y[i] - x[i]) / tol); 3559 } 3560 } 3561 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 3562 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 3563 3564 ierr = MPI_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,((PetscObject)ts)->comm);CHKERRQ(ierr); 3565 *norm = PetscSqrtReal(gsum / N); 3566 if (PetscIsInfOrNanScalar(*norm)) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 3567 PetscFunctionReturn(0); 3568 } 3569 3570 #undef __FUNCT__ 3571 #define __FUNCT__ "TSSetCFLTimeLocal" 3572 /*@ 3573 TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler 3574 3575 Logically Collective on TS 3576 3577 Input Arguments: 3578 + ts - time stepping context 3579 - cfltime - maximum stable time step if using forward Euler (value can be different on each process) 3580 3581 Note: 3582 After calling this function, the global CFL time can be obtained by calling TSGetCFLTime() 3583 3584 Level: intermediate 3585 3586 .seealso: TSGetCFLTime(), TSADAPTCFL 3587 @*/ 3588 PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime) 3589 { 3590 3591 PetscFunctionBegin; 3592 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3593 ts->cfltime_local = cfltime; 3594 ts->cfltime = -1.; 3595 PetscFunctionReturn(0); 3596 } 3597 3598 #undef __FUNCT__ 3599 #define __FUNCT__ "TSGetCFLTime" 3600 /*@ 3601 TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler 3602 3603 Collective on TS 3604 3605 Input Arguments: 3606 . ts - time stepping context 3607 3608 Output Arguments: 3609 . cfltime - maximum stable time step for forward Euler 3610 3611 Level: advanced 3612 3613 .seealso: TSSetCFLTimeLocal() 3614 @*/ 3615 PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime) 3616 { 3617 PetscErrorCode ierr; 3618 3619 PetscFunctionBegin; 3620 if (ts->cfltime < 0) { 3621 ierr = MPI_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,((PetscObject)ts)->comm);CHKERRQ(ierr); 3622 } 3623 *cfltime = ts->cfltime; 3624 PetscFunctionReturn(0); 3625 } 3626 3627 #undef __FUNCT__ 3628 #define __FUNCT__ "TSVISetVariableBounds" 3629 /*@ 3630 TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu 3631 3632 Input Parameters: 3633 . ts - the TS context. 3634 . xl - lower bound. 3635 . xu - upper bound. 3636 3637 Notes: 3638 If this routine is not called then the lower and upper bounds are set to 3639 SNES_VI_NINF and SNES_VI_INF respectively during SNESSetUp(). 3640 3641 Level: advanced 3642 3643 @*/ 3644 PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu) 3645 { 3646 PetscErrorCode ierr; 3647 SNES snes; 3648 3649 PetscFunctionBegin; 3650 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3651 ierr = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(ierr); 3652 PetscFunctionReturn(0); 3653 } 3654 3655 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3656 #include <mex.h> 3657 3658 typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext; 3659 3660 #undef __FUNCT__ 3661 #define __FUNCT__ "TSComputeFunction_Matlab" 3662 /* 3663 TSComputeFunction_Matlab - Calls the function that has been set with 3664 TSSetFunctionMatlab(). 3665 3666 Collective on TS 3667 3668 Input Parameters: 3669 + snes - the TS context 3670 - x - input vector 3671 3672 Output Parameter: 3673 . y - function vector, as set by TSSetFunction() 3674 3675 Notes: 3676 TSComputeFunction() is typically used within nonlinear solvers 3677 implementations, so most users would not generally call this routine 3678 themselves. 3679 3680 Level: developer 3681 3682 .keywords: TS, nonlinear, compute, function 3683 3684 .seealso: TSSetFunction(), TSGetFunction() 3685 */ 3686 PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec x,Vec xdot,Vec y, void *ctx) 3687 { 3688 PetscErrorCode ierr; 3689 TSMatlabContext *sctx = (TSMatlabContext *)ctx; 3690 int nlhs = 1,nrhs = 7; 3691 mxArray *plhs[1],*prhs[7]; 3692 long long int lx = 0,lxdot = 0,ly = 0,ls = 0; 3693 3694 PetscFunctionBegin; 3695 PetscValidHeaderSpecific(snes,TS_CLASSID,1); 3696 PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3697 PetscValidHeaderSpecific(xdot,VEC_CLASSID,4); 3698 PetscValidHeaderSpecific(y,VEC_CLASSID,5); 3699 PetscCheckSameComm(snes,1,x,3); 3700 PetscCheckSameComm(snes,1,y,5); 3701 3702 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 3703 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 3704 ierr = PetscMemcpy(&lxdot,&xdot,sizeof(xdot));CHKERRQ(ierr); 3705 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 3706 prhs[0] = mxCreateDoubleScalar((double)ls); 3707 prhs[1] = mxCreateDoubleScalar(time); 3708 prhs[2] = mxCreateDoubleScalar((double)lx); 3709 prhs[3] = mxCreateDoubleScalar((double)lxdot); 3710 prhs[4] = mxCreateDoubleScalar((double)ly); 3711 prhs[5] = mxCreateString(sctx->funcname); 3712 prhs[6] = sctx->ctx; 3713 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");CHKERRQ(ierr); 3714 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 3715 mxDestroyArray(prhs[0]); 3716 mxDestroyArray(prhs[1]); 3717 mxDestroyArray(prhs[2]); 3718 mxDestroyArray(prhs[3]); 3719 mxDestroyArray(prhs[4]); 3720 mxDestroyArray(prhs[5]); 3721 mxDestroyArray(plhs[0]); 3722 PetscFunctionReturn(0); 3723 } 3724 3725 3726 #undef __FUNCT__ 3727 #define __FUNCT__ "TSSetFunctionMatlab" 3728 /* 3729 TSSetFunctionMatlab - Sets the function evaluation routine and function 3730 vector for use by the TS routines in solving ODEs 3731 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 3732 3733 Logically Collective on TS 3734 3735 Input Parameters: 3736 + ts - the TS context 3737 - func - function evaluation routine 3738 3739 Calling sequence of func: 3740 $ func (TS ts,PetscReal time,Vec x,Vec xdot,Vec f,void *ctx); 3741 3742 Level: beginner 3743 3744 .keywords: TS, nonlinear, set, function 3745 3746 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 3747 */ 3748 PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx) 3749 { 3750 PetscErrorCode ierr; 3751 TSMatlabContext *sctx; 3752 3753 PetscFunctionBegin; 3754 /* currently sctx is memory bleed */ 3755 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 3756 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 3757 /* 3758 This should work, but it doesn't 3759 sctx->ctx = ctx; 3760 mexMakeArrayPersistent(sctx->ctx); 3761 */ 3762 sctx->ctx = mxDuplicateArray(ctx); 3763 ierr = TSSetIFunction(ts,PETSC_NULL,TSComputeFunction_Matlab,sctx);CHKERRQ(ierr); 3764 PetscFunctionReturn(0); 3765 } 3766 3767 #undef __FUNCT__ 3768 #define __FUNCT__ "TSComputeJacobian_Matlab" 3769 /* 3770 TSComputeJacobian_Matlab - Calls the function that has been set with 3771 TSSetJacobianMatlab(). 3772 3773 Collective on TS 3774 3775 Input Parameters: 3776 + ts - the TS context 3777 . x - input vector 3778 . A, B - the matrices 3779 - ctx - user context 3780 3781 Output Parameter: 3782 . flag - structure of the matrix 3783 3784 Level: developer 3785 3786 .keywords: TS, nonlinear, compute, function 3787 3788 .seealso: TSSetFunction(), TSGetFunction() 3789 @*/ 3790 PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec x,Vec xdot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag, void *ctx) 3791 { 3792 PetscErrorCode ierr; 3793 TSMatlabContext *sctx = (TSMatlabContext *)ctx; 3794 int nlhs = 2,nrhs = 9; 3795 mxArray *plhs[2],*prhs[9]; 3796 long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0; 3797 3798 PetscFunctionBegin; 3799 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3800 PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3801 3802 /* call Matlab function in ctx with arguments x and y */ 3803 3804 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 3805 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 3806 ierr = PetscMemcpy(&lxdot,&xdot,sizeof(x));CHKERRQ(ierr); 3807 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 3808 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 3809 prhs[0] = mxCreateDoubleScalar((double)ls); 3810 prhs[1] = mxCreateDoubleScalar((double)time); 3811 prhs[2] = mxCreateDoubleScalar((double)lx); 3812 prhs[3] = mxCreateDoubleScalar((double)lxdot); 3813 prhs[4] = mxCreateDoubleScalar((double)shift); 3814 prhs[5] = mxCreateDoubleScalar((double)lA); 3815 prhs[6] = mxCreateDoubleScalar((double)lB); 3816 prhs[7] = mxCreateString(sctx->funcname); 3817 prhs[8] = sctx->ctx; 3818 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");CHKERRQ(ierr); 3819 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 3820 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 3821 mxDestroyArray(prhs[0]); 3822 mxDestroyArray(prhs[1]); 3823 mxDestroyArray(prhs[2]); 3824 mxDestroyArray(prhs[3]); 3825 mxDestroyArray(prhs[4]); 3826 mxDestroyArray(prhs[5]); 3827 mxDestroyArray(prhs[6]); 3828 mxDestroyArray(prhs[7]); 3829 mxDestroyArray(plhs[0]); 3830 mxDestroyArray(plhs[1]); 3831 PetscFunctionReturn(0); 3832 } 3833 3834 3835 #undef __FUNCT__ 3836 #define __FUNCT__ "TSSetJacobianMatlab" 3837 /* 3838 TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 3839 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 3840 3841 Logically Collective on TS 3842 3843 Input Parameters: 3844 + ts - the TS context 3845 . A,B - Jacobian matrices 3846 . func - function evaluation routine 3847 - ctx - user context 3848 3849 Calling sequence of func: 3850 $ flag = func (TS ts,PetscReal time,Vec x,Vec xdot,Mat A,Mat B,void *ctx); 3851 3852 3853 Level: developer 3854 3855 .keywords: TS, nonlinear, set, function 3856 3857 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 3858 */ 3859 PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx) 3860 { 3861 PetscErrorCode ierr; 3862 TSMatlabContext *sctx; 3863 3864 PetscFunctionBegin; 3865 /* currently sctx is memory bleed */ 3866 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 3867 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 3868 /* 3869 This should work, but it doesn't 3870 sctx->ctx = ctx; 3871 mexMakeArrayPersistent(sctx->ctx); 3872 */ 3873 sctx->ctx = mxDuplicateArray(ctx); 3874 ierr = TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 3875 PetscFunctionReturn(0); 3876 } 3877 3878 #undef __FUNCT__ 3879 #define __FUNCT__ "TSMonitor_Matlab" 3880 /* 3881 TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab(). 3882 3883 Collective on TS 3884 3885 .seealso: TSSetFunction(), TSGetFunction() 3886 @*/ 3887 PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec x, void *ctx) 3888 { 3889 PetscErrorCode ierr; 3890 TSMatlabContext *sctx = (TSMatlabContext *)ctx; 3891 int nlhs = 1,nrhs = 6; 3892 mxArray *plhs[1],*prhs[6]; 3893 long long int lx = 0,ls = 0; 3894 3895 PetscFunctionBegin; 3896 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3897 PetscValidHeaderSpecific(x,VEC_CLASSID,4); 3898 3899 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 3900 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 3901 prhs[0] = mxCreateDoubleScalar((double)ls); 3902 prhs[1] = mxCreateDoubleScalar((double)it); 3903 prhs[2] = mxCreateDoubleScalar((double)time); 3904 prhs[3] = mxCreateDoubleScalar((double)lx); 3905 prhs[4] = mxCreateString(sctx->funcname); 3906 prhs[5] = sctx->ctx; 3907 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");CHKERRQ(ierr); 3908 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 3909 mxDestroyArray(prhs[0]); 3910 mxDestroyArray(prhs[1]); 3911 mxDestroyArray(prhs[2]); 3912 mxDestroyArray(prhs[3]); 3913 mxDestroyArray(prhs[4]); 3914 mxDestroyArray(plhs[0]); 3915 PetscFunctionReturn(0); 3916 } 3917 3918 3919 #undef __FUNCT__ 3920 #define __FUNCT__ "TSMonitorSetMatlab" 3921 /* 3922 TSMonitorSetMatlab - Sets the monitor function from Matlab 3923 3924 Level: developer 3925 3926 .keywords: TS, nonlinear, set, function 3927 3928 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 3929 */ 3930 PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx) 3931 { 3932 PetscErrorCode ierr; 3933 TSMatlabContext *sctx; 3934 3935 PetscFunctionBegin; 3936 /* currently sctx is memory bleed */ 3937 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 3938 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 3939 /* 3940 This should work, but it doesn't 3941 sctx->ctx = ctx; 3942 mexMakeArrayPersistent(sctx->ctx); 3943 */ 3944 sctx->ctx = mxDuplicateArray(ctx); 3945 ierr = TSMonitorSet(ts,TSMonitor_Matlab,sctx,PETSC_NULL);CHKERRQ(ierr); 3946 PetscFunctionReturn(0); 3947 } 3948 #endif 3949 3950 #undef __FUNCT__ 3951 #define __FUNCT__ "TSMonitorSolutionODE" 3952 /*@C 3953 TSMonitorSolutionODE - Monitors progress of the TS solvers by plotting each component of the solution vector 3954 in a time based line graph 3955 3956 Collective on TS 3957 3958 Input Parameters: 3959 + ts - the TS context 3960 . step - current time-step 3961 . ptime - current time 3962 - lg - a line graph object 3963 3964 Level: intermediate 3965 3966 Notes: only for sequential solves 3967 3968 .keywords: TS, vector, monitor, view 3969 3970 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3971 @*/ 3972 PetscErrorCode TSMonitorSolutionODE(TS ts,PetscInt step,PetscReal ptime,Vec x,void *dummy) 3973 { 3974 PetscErrorCode ierr; 3975 PetscDrawLG lg = (PetscDrawLG)dummy; 3976 const PetscScalar *yy; 3977 3978 PetscFunctionBegin; 3979 if (!step) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3980 ierr = VecGetArrayRead(x,&yy);CHKERRQ(ierr); 3981 ierr = PetscDrawLGAddCommonPoint(lg,ptime,yy);CHKERRQ(ierr); 3982 ierr = VecRestoreArrayRead(x,&yy);CHKERRQ(ierr); 3983 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3984 PetscFunctionReturn(0); 3985 } 3986 3987 #undef __FUNCT__ 3988 #define __FUNCT__ "TSMonitorSolutionODEDestroy" 3989 /*@C 3990 TSMonitorSolutionODEDestroy - Destroys the monitor context for TSMonitorSolutionODE() 3991 3992 Collective on TS 3993 3994 Input Parameters: 3995 . ctx - the monitor context 3996 3997 Level: intermediate 3998 3999 .keywords: TS, vector, monitor, view 4000 4001 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorSolutionODE() 4002 @*/ 4003 PetscErrorCode TSMonitorSolutionODEDestroy(PetscDrawLG *lg) 4004 { 4005 PetscDraw draw; 4006 PetscErrorCode ierr; 4007 4008 PetscFunctionBegin; 4009 ierr = PetscDrawLGGetDraw(*lg,&draw);CHKERRQ(ierr); 4010 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 4011 ierr = PetscDrawLGDestroy(lg);CHKERRQ(ierr); 4012 PetscFunctionReturn(0); 4013 } 4014 4015 #undef __FUNCT__ 4016 #define __FUNCT__ "TSMonitorSolutionODECreate" 4017 /*@C 4018 TSMonitorSolutionODECreate - Creates the monitor context for TSMonitorSolutionODE() 4019 4020 Collective on TS 4021 4022 Input Parameter: 4023 + comm - MPI communicator 4024 . N - number of components in the solution vector 4025 - 4026 4027 Output Patameter: 4028 . ctx - the monitor context 4029 4030 Level: intermediate 4031 4032 .keywords: TS, vector, monitor, view 4033 4034 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorSolution() 4035 @*/ 4036 PetscErrorCode TSMonitorSolutionODECreate(MPI_Comm comm,PetscInt N,const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 4037 { 4038 PetscDraw win; 4039 PetscErrorCode ierr; 4040 PetscDrawAxis axis; 4041 4042 PetscFunctionBegin; 4043 ierr = PetscDrawCreate(comm,host,label,x,y,m,n,&win);CHKERRQ(ierr); 4044 ierr = PetscDrawSetType(win,PETSC_DRAW_X);CHKERRQ(ierr); 4045 ierr = PetscDrawLGCreate(win,N,draw);CHKERRQ(ierr); 4046 ierr = PetscDrawLGIndicateDataPoints(*draw);CHKERRQ(ierr); 4047 ierr = PetscDrawLGGetAxis(*draw,&axis);CHKERRQ(ierr); 4048 ierr = PetscDrawAxisSetLabels(axis,"Solution","Time","Solution");CHKERRQ(ierr); 4049 ierr = PetscLogObjectParent(*draw,win);CHKERRQ(ierr); 4050 PetscFunctionReturn(0); 4051 } 4052 4053 #undef __FUNCT__ 4054 #define __FUNCT__ "TSMonitorErrorODE" 4055 /*@C 4056 TSMonitorErrorODE - Monitors progress of the TS solvers by plotting each component of the solution vector 4057 in a time based line graph 4058 4059 Collective on TS 4060 4061 Input Parameters: 4062 + ts - the TS context 4063 . step - current time-step 4064 . ptime - current time 4065 - lg - a line graph object 4066 4067 Level: intermediate 4068 4069 Notes: 4070 Only for sequential solves. 4071 4072 The user must provide the solution using TSSetSolutionFunction() to use this monitor. 4073 4074 Options Database Keys: 4075 . -ts_monitor_errorode - create a graphical monitor of error history 4076 4077 .keywords: TS, vector, monitor, view 4078 4079 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction() 4080 @*/ 4081 PetscErrorCode TSMonitorErrorODE(TS ts,PetscInt step,PetscReal ptime,Vec x,void *dummy) 4082 { 4083 PetscErrorCode ierr; 4084 PetscDrawLG lg = (PetscDrawLG)dummy; 4085 const PetscScalar *yy; 4086 Vec y; 4087 4088 PetscFunctionBegin; 4089 if (!step) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 4090 ierr = VecDuplicate(x,&y);CHKERRQ(ierr); 4091 ierr = TSComputeSolutionFunction(ts,ptime,y);CHKERRQ(ierr); 4092 ierr = VecAXPY(y,-1.0,x);CHKERRQ(ierr); 4093 ierr = VecGetArrayRead(y,&yy);CHKERRQ(ierr); 4094 ierr = PetscDrawLGAddCommonPoint(lg,ptime,yy);CHKERRQ(ierr); 4095 ierr = VecRestoreArrayRead(y,&yy);CHKERRQ(ierr); 4096 ierr = VecDestroy(&y);CHKERRQ(ierr); 4097 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 4098 PetscFunctionReturn(0); 4099 } 4100 4101 #undef __FUNCT__ 4102 #define __FUNCT__ "TSMonitorErrorODEDestroy" 4103 /*@C 4104 TSMonitorErrorODEDestroy - Destroys the monitor context for TSMonitorErrorODE() 4105 4106 Collective on TS 4107 4108 Input Parameters: 4109 . ctx - the monitor context 4110 4111 Level: intermediate 4112 4113 .keywords: TS, vector, monitor, view 4114 4115 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorErrorODE() 4116 @*/ 4117 PetscErrorCode TSMonitorErrorODEDestroy(PetscDrawLG *lg) 4118 { 4119 PetscDraw draw; 4120 PetscErrorCode ierr; 4121 4122 PetscFunctionBegin; 4123 ierr = PetscDrawLGGetDraw(*lg,&draw);CHKERRQ(ierr); 4124 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 4125 ierr = PetscDrawLGDestroy(lg);CHKERRQ(ierr); 4126 PetscFunctionReturn(0); 4127 } 4128 4129 #undef __FUNCT__ 4130 #define __FUNCT__ "TSMonitorErrorODECreate" 4131 /*@C 4132 TSMonitorErrorODECreate - Creates the monitor context for TSMonitorErrorODE() 4133 4134 Collective on TS 4135 4136 Input Parameter: 4137 + comm - MPI communicator 4138 . N - number of components in the solution vector 4139 - 4140 4141 Output Patameter: 4142 . ctx - the monitor context 4143 4144 Options Database Keys: 4145 . -ts_monitor_errorode - create a graphical monitor of error history 4146 4147 Level: intermediate 4148 4149 .keywords: TS, vector, monitor, view 4150 4151 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorErrorODE() 4152 @*/ 4153 PetscErrorCode TSMonitorErrorODECreate(MPI_Comm comm,PetscInt N,const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 4154 { 4155 PetscDraw win; 4156 PetscErrorCode ierr; 4157 PetscDrawAxis axis; 4158 4159 PetscFunctionBegin; 4160 ierr = PetscDrawCreate(comm,host,label,x,y,m,n,&win);CHKERRQ(ierr); 4161 ierr = PetscDrawSetType(win,PETSC_DRAW_X);CHKERRQ(ierr); 4162 ierr = PetscDrawLGCreate(win,N,draw);CHKERRQ(ierr); 4163 ierr = PetscDrawLGIndicateDataPoints(*draw);CHKERRQ(ierr); 4164 ierr = PetscDrawLGGetAxis(*draw,&axis);CHKERRQ(ierr); 4165 ierr = PetscDrawAxisSetLabels(axis,"Error","Time","Error");CHKERRQ(ierr); 4166 ierr = PetscLogObjectParent(*draw,win);CHKERRQ(ierr); 4167 PetscFunctionReturn(0); 4168 } 4169