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