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","TSSetMaxSNESFailures",ts->max_snes_failures,&ts->max_snes_failures,PETSC_NULL);CHKERRQ(ierr); 96 ierr = PetscOptionsInt("-ts_max_reject","Maximum number of step rejections before step fails","TSSetMaxStepRejections",ts->max_reject,&ts->max_reject,PETSC_NULL);CHKERRQ(ierr); 97 ierr = PetscOptionsBool("-ts_error_if_step_fails","Error if no step succeeds","TSSetErrorIfStepFails",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->snes_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->ksp_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 = VecDestroy(&ts->vatol);CHKERRQ(ierr); 1267 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 1268 ierr = VecDestroyVecs(ts->nwork,&ts->work);CHKERRQ(ierr); 1269 ts->setupcalled = PETSC_FALSE; 1270 PetscFunctionReturn(0); 1271 } 1272 1273 #undef __FUNCT__ 1274 #define __FUNCT__ "TSDestroy" 1275 /*@ 1276 TSDestroy - Destroys the timestepper context that was created 1277 with TSCreate(). 1278 1279 Collective on TS 1280 1281 Input Parameter: 1282 . ts - the TS context obtained from TSCreate() 1283 1284 Level: beginner 1285 1286 .keywords: TS, timestepper, destroy 1287 1288 .seealso: TSCreate(), TSSetUp(), TSSolve() 1289 @*/ 1290 PetscErrorCode TSDestroy(TS *ts) 1291 { 1292 PetscErrorCode ierr; 1293 1294 PetscFunctionBegin; 1295 if (!*ts) PetscFunctionReturn(0); 1296 PetscValidHeaderSpecific((*ts),TS_CLASSID,1); 1297 if (--((PetscObject)(*ts))->refct > 0) {*ts = 0; PetscFunctionReturn(0);} 1298 1299 ierr = TSReset((*ts));CHKERRQ(ierr); 1300 1301 /* if memory was published with AMS then destroy it */ 1302 ierr = PetscObjectDepublish((*ts));CHKERRQ(ierr); 1303 if ((*ts)->ops->destroy) {ierr = (*(*ts)->ops->destroy)((*ts));CHKERRQ(ierr);} 1304 1305 ierr = TSAdaptDestroy(&(*ts)->adapt);CHKERRQ(ierr); 1306 ierr = SNESDestroy(&(*ts)->snes);CHKERRQ(ierr); 1307 ierr = DMDestroy(&(*ts)->dm);CHKERRQ(ierr); 1308 ierr = TSMonitorCancel((*ts));CHKERRQ(ierr); 1309 1310 ierr = PetscFree((*ts)->userops); 1311 1312 ierr = PetscHeaderDestroy(ts);CHKERRQ(ierr); 1313 PetscFunctionReturn(0); 1314 } 1315 1316 #undef __FUNCT__ 1317 #define __FUNCT__ "TSGetSNES" 1318 /*@ 1319 TSGetSNES - Returns the SNES (nonlinear solver) associated with 1320 a TS (timestepper) context. Valid only for nonlinear problems. 1321 1322 Not Collective, but SNES is parallel if TS is parallel 1323 1324 Input Parameter: 1325 . ts - the TS context obtained from TSCreate() 1326 1327 Output Parameter: 1328 . snes - the nonlinear solver context 1329 1330 Notes: 1331 The user can then directly manipulate the SNES context to set various 1332 options, etc. Likewise, the user can then extract and manipulate the 1333 KSP, KSP, and PC contexts as well. 1334 1335 TSGetSNES() does not work for integrators that do not use SNES; in 1336 this case TSGetSNES() returns PETSC_NULL in snes. 1337 1338 Level: beginner 1339 1340 .keywords: timestep, get, SNES 1341 @*/ 1342 PetscErrorCode TSGetSNES(TS ts,SNES *snes) 1343 { 1344 PetscErrorCode ierr; 1345 1346 PetscFunctionBegin; 1347 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1348 PetscValidPointer(snes,2); 1349 if (!ts->snes) { 1350 ierr = SNESCreate(((PetscObject)ts)->comm,&ts->snes);CHKERRQ(ierr); 1351 ierr = PetscLogObjectParent(ts,ts->snes);CHKERRQ(ierr); 1352 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);CHKERRQ(ierr); 1353 if (ts->dm) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 1354 if (ts->problem_type == TS_LINEAR) { 1355 ierr = SNESSetType(ts->snes,SNESKSPONLY);CHKERRQ(ierr); 1356 } 1357 } 1358 *snes = ts->snes; 1359 PetscFunctionReturn(0); 1360 } 1361 1362 #undef __FUNCT__ 1363 #define __FUNCT__ "TSGetKSP" 1364 /*@ 1365 TSGetKSP - Returns the KSP (linear solver) associated with 1366 a TS (timestepper) context. 1367 1368 Not Collective, but KSP is parallel if TS is parallel 1369 1370 Input Parameter: 1371 . ts - the TS context obtained from TSCreate() 1372 1373 Output Parameter: 1374 . ksp - the nonlinear solver context 1375 1376 Notes: 1377 The user can then directly manipulate the KSP context to set various 1378 options, etc. Likewise, the user can then extract and manipulate the 1379 KSP and PC contexts as well. 1380 1381 TSGetKSP() does not work for integrators that do not use KSP; 1382 in this case TSGetKSP() returns PETSC_NULL in ksp. 1383 1384 Level: beginner 1385 1386 .keywords: timestep, get, KSP 1387 @*/ 1388 PetscErrorCode TSGetKSP(TS ts,KSP *ksp) 1389 { 1390 PetscErrorCode ierr; 1391 SNES snes; 1392 1393 PetscFunctionBegin; 1394 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1395 PetscValidPointer(ksp,2); 1396 if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first"); 1397 if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()"); 1398 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1399 ierr = SNESGetKSP(snes,ksp);CHKERRQ(ierr); 1400 PetscFunctionReturn(0); 1401 } 1402 1403 /* ----------- Routines to set solver parameters ---------- */ 1404 1405 #undef __FUNCT__ 1406 #define __FUNCT__ "TSGetDuration" 1407 /*@ 1408 TSGetDuration - Gets the maximum number of timesteps to use and 1409 maximum time for iteration. 1410 1411 Not Collective 1412 1413 Input Parameters: 1414 + ts - the TS context obtained from TSCreate() 1415 . maxsteps - maximum number of iterations to use, or PETSC_NULL 1416 - maxtime - final time to iterate to, or PETSC_NULL 1417 1418 Level: intermediate 1419 1420 .keywords: TS, timestep, get, maximum, iterations, time 1421 @*/ 1422 PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime) 1423 { 1424 PetscFunctionBegin; 1425 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1426 if (maxsteps) { 1427 PetscValidIntPointer(maxsteps,2); 1428 *maxsteps = ts->max_steps; 1429 } 1430 if (maxtime) { 1431 PetscValidScalarPointer(maxtime,3); 1432 *maxtime = ts->max_time; 1433 } 1434 PetscFunctionReturn(0); 1435 } 1436 1437 #undef __FUNCT__ 1438 #define __FUNCT__ "TSSetDuration" 1439 /*@ 1440 TSSetDuration - Sets the maximum number of timesteps to use and 1441 maximum time for iteration. 1442 1443 Logically Collective on TS 1444 1445 Input Parameters: 1446 + ts - the TS context obtained from TSCreate() 1447 . maxsteps - maximum number of iterations to use 1448 - maxtime - final time to iterate to 1449 1450 Options Database Keys: 1451 . -ts_max_steps <maxsteps> - Sets maxsteps 1452 . -ts_final_time <maxtime> - Sets maxtime 1453 1454 Notes: 1455 The default maximum number of iterations is 5000. Default time is 5.0 1456 1457 Level: intermediate 1458 1459 .keywords: TS, timestep, set, maximum, iterations 1460 1461 .seealso: TSSetExactFinalTime() 1462 @*/ 1463 PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime) 1464 { 1465 PetscFunctionBegin; 1466 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1467 PetscValidLogicalCollectiveInt(ts,maxsteps,2); 1468 PetscValidLogicalCollectiveReal(ts,maxtime,2); 1469 if (maxsteps >= 0) ts->max_steps = maxsteps; 1470 if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime; 1471 PetscFunctionReturn(0); 1472 } 1473 1474 #undef __FUNCT__ 1475 #define __FUNCT__ "TSSetSolution" 1476 /*@ 1477 TSSetSolution - Sets the initial solution vector 1478 for use by the TS routines. 1479 1480 Logically Collective on TS and Vec 1481 1482 Input Parameters: 1483 + ts - the TS context obtained from TSCreate() 1484 - x - the solution vector 1485 1486 Level: beginner 1487 1488 .keywords: TS, timestep, set, solution, initial conditions 1489 @*/ 1490 PetscErrorCode TSSetSolution(TS ts,Vec x) 1491 { 1492 PetscErrorCode ierr; 1493 DM dm; 1494 1495 PetscFunctionBegin; 1496 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1497 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 1498 ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr); 1499 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 1500 ts->vec_sol = x; 1501 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1502 ierr = DMShellSetGlobalVector(dm,x);CHKERRQ(ierr); 1503 PetscFunctionReturn(0); 1504 } 1505 1506 #undef __FUNCT__ 1507 #define __FUNCT__ "TSSetPreStep" 1508 /*@C 1509 TSSetPreStep - Sets the general-purpose function 1510 called once at the beginning of each time step. 1511 1512 Logically Collective on TS 1513 1514 Input Parameters: 1515 + ts - The TS context obtained from TSCreate() 1516 - func - The function 1517 1518 Calling sequence of func: 1519 . func (TS ts); 1520 1521 Level: intermediate 1522 1523 .keywords: TS, timestep 1524 @*/ 1525 PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS)) 1526 { 1527 PetscFunctionBegin; 1528 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1529 ts->ops->prestep = func; 1530 PetscFunctionReturn(0); 1531 } 1532 1533 #undef __FUNCT__ 1534 #define __FUNCT__ "TSPreStep" 1535 /*@ 1536 TSPreStep - Runs the user-defined pre-step function. 1537 1538 Collective on TS 1539 1540 Input Parameters: 1541 . ts - The TS context obtained from TSCreate() 1542 1543 Notes: 1544 TSPreStep() is typically used within time stepping implementations, 1545 so most users would not generally call this routine themselves. 1546 1547 Level: developer 1548 1549 .keywords: TS, timestep 1550 @*/ 1551 PetscErrorCode TSPreStep(TS ts) 1552 { 1553 PetscErrorCode ierr; 1554 1555 PetscFunctionBegin; 1556 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1557 if (ts->ops->prestep) { 1558 PetscStackPush("TS PreStep function"); 1559 ierr = (*ts->ops->prestep)(ts);CHKERRQ(ierr); 1560 PetscStackPop; 1561 } 1562 PetscFunctionReturn(0); 1563 } 1564 1565 #undef __FUNCT__ 1566 #define __FUNCT__ "TSSetPostStep" 1567 /*@C 1568 TSSetPostStep - Sets the general-purpose function 1569 called once at the end of each time step. 1570 1571 Logically Collective on TS 1572 1573 Input Parameters: 1574 + ts - The TS context obtained from TSCreate() 1575 - func - The function 1576 1577 Calling sequence of func: 1578 . func (TS ts); 1579 1580 Level: intermediate 1581 1582 .keywords: TS, timestep 1583 @*/ 1584 PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS)) 1585 { 1586 PetscFunctionBegin; 1587 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1588 ts->ops->poststep = func; 1589 PetscFunctionReturn(0); 1590 } 1591 1592 #undef __FUNCT__ 1593 #define __FUNCT__ "TSPostStep" 1594 /*@ 1595 TSPostStep - Runs the user-defined post-step function. 1596 1597 Collective on TS 1598 1599 Input Parameters: 1600 . ts - The TS context obtained from TSCreate() 1601 1602 Notes: 1603 TSPostStep() is typically used within time stepping implementations, 1604 so most users would not generally call this routine themselves. 1605 1606 Level: developer 1607 1608 .keywords: TS, timestep 1609 @*/ 1610 PetscErrorCode TSPostStep(TS ts) 1611 { 1612 PetscErrorCode ierr; 1613 1614 PetscFunctionBegin; 1615 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1616 if (ts->ops->poststep) { 1617 PetscStackPush("TS PostStep function"); 1618 ierr = (*ts->ops->poststep)(ts);CHKERRQ(ierr); 1619 PetscStackPop; 1620 } 1621 PetscFunctionReturn(0); 1622 } 1623 1624 /* ------------ Routines to set performance monitoring options ----------- */ 1625 1626 #undef __FUNCT__ 1627 #define __FUNCT__ "TSMonitorSet" 1628 /*@C 1629 TSMonitorSet - Sets an ADDITIONAL function that is to be used at every 1630 timestep to display the iteration's progress. 1631 1632 Logically Collective on TS 1633 1634 Input Parameters: 1635 + ts - the TS context obtained from TSCreate() 1636 . monitor - monitoring routine 1637 . mctx - [optional] user-defined context for private data for the 1638 monitor routine (use PETSC_NULL if no context is desired) 1639 - monitordestroy - [optional] routine that frees monitor context 1640 (may be PETSC_NULL) 1641 1642 Calling sequence of monitor: 1643 $ int monitor(TS ts,PetscInt steps,PetscReal time,Vec x,void *mctx) 1644 1645 + ts - the TS context 1646 . steps - iteration number 1647 . time - current time 1648 . x - current iterate 1649 - mctx - [optional] monitoring context 1650 1651 Notes: 1652 This routine adds an additional monitor to the list of monitors that 1653 already has been loaded. 1654 1655 Fortran notes: Only a single monitor function can be set for each TS object 1656 1657 Level: intermediate 1658 1659 .keywords: TS, timestep, set, monitor 1660 1661 .seealso: TSMonitorDefault(), TSMonitorCancel() 1662 @*/ 1663 PetscErrorCode TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**)) 1664 { 1665 PetscFunctionBegin; 1666 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1667 if (ts->numbermonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 1668 ts->monitor[ts->numbermonitors] = monitor; 1669 ts->mdestroy[ts->numbermonitors] = mdestroy; 1670 ts->monitorcontext[ts->numbermonitors++] = (void*)mctx; 1671 PetscFunctionReturn(0); 1672 } 1673 1674 #undef __FUNCT__ 1675 #define __FUNCT__ "TSMonitorCancel" 1676 /*@C 1677 TSMonitorCancel - Clears all the monitors that have been set on a time-step object. 1678 1679 Logically Collective on TS 1680 1681 Input Parameters: 1682 . ts - the TS context obtained from TSCreate() 1683 1684 Notes: 1685 There is no way to remove a single, specific monitor. 1686 1687 Level: intermediate 1688 1689 .keywords: TS, timestep, set, monitor 1690 1691 .seealso: TSMonitorDefault(), TSMonitorSet() 1692 @*/ 1693 PetscErrorCode TSMonitorCancel(TS ts) 1694 { 1695 PetscErrorCode ierr; 1696 PetscInt i; 1697 1698 PetscFunctionBegin; 1699 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1700 for (i=0; i<ts->numbermonitors; i++) { 1701 if (ts->mdestroy[i]) { 1702 ierr = (*ts->mdestroy[i])(&ts->monitorcontext[i]);CHKERRQ(ierr); 1703 } 1704 } 1705 ts->numbermonitors = 0; 1706 PetscFunctionReturn(0); 1707 } 1708 1709 #undef __FUNCT__ 1710 #define __FUNCT__ "TSMonitorDefault" 1711 /*@ 1712 TSMonitorDefault - Sets the Default monitor 1713 1714 Level: intermediate 1715 1716 .keywords: TS, set, monitor 1717 1718 .seealso: TSMonitorDefault(), TSMonitorSet() 1719 @*/ 1720 PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy) 1721 { 1722 PetscErrorCode ierr; 1723 PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(((PetscObject)ts)->comm); 1724 1725 PetscFunctionBegin; 1726 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 1727 ierr = PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g\n",step,(double)ts->time_step,(double)ptime);CHKERRQ(ierr); 1728 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 1729 PetscFunctionReturn(0); 1730 } 1731 1732 #undef __FUNCT__ 1733 #define __FUNCT__ "TSSetRetainStages" 1734 /*@ 1735 TSSetRetainStages - Request that all stages in the upcoming step be stored so that interpolation will be available. 1736 1737 Logically Collective on TS 1738 1739 Input Argument: 1740 . ts - time stepping context 1741 1742 Output Argument: 1743 . flg - PETSC_TRUE or PETSC_FALSE 1744 1745 Level: intermediate 1746 1747 .keywords: TS, set 1748 1749 .seealso: TSInterpolate(), TSSetPostStep() 1750 @*/ 1751 PetscErrorCode TSSetRetainStages(TS ts,PetscBool flg) 1752 { 1753 1754 PetscFunctionBegin; 1755 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1756 ts->retain_stages = flg; 1757 PetscFunctionReturn(0); 1758 } 1759 1760 #undef __FUNCT__ 1761 #define __FUNCT__ "TSInterpolate" 1762 /*@ 1763 TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval 1764 1765 Collective on TS 1766 1767 Input Argument: 1768 + ts - time stepping context 1769 - t - time to interpolate to 1770 1771 Output Argument: 1772 . X - state at given time 1773 1774 Notes: 1775 The user should call TSSetRetainStages() before taking a step in which interpolation will be requested. 1776 1777 Level: intermediate 1778 1779 Developer Notes: 1780 TSInterpolate() and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints. 1781 1782 .keywords: TS, set 1783 1784 .seealso: TSSetRetainStages(), TSSetPostStep() 1785 @*/ 1786 PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec X) 1787 { 1788 PetscErrorCode ierr; 1789 1790 PetscFunctionBegin; 1791 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1792 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); 1793 if (!ts->ops->interpolate) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name); 1794 ierr = (*ts->ops->interpolate)(ts,t,X);CHKERRQ(ierr); 1795 PetscFunctionReturn(0); 1796 } 1797 1798 #undef __FUNCT__ 1799 #define __FUNCT__ "TSStep" 1800 /*@ 1801 TSStep - Steps one time step 1802 1803 Collective on TS 1804 1805 Input Parameter: 1806 . ts - the TS context obtained from TSCreate() 1807 1808 Level: intermediate 1809 1810 .keywords: TS, timestep, solve 1811 1812 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve() 1813 @*/ 1814 PetscErrorCode TSStep(TS ts) 1815 { 1816 PetscReal ptime_prev; 1817 PetscErrorCode ierr; 1818 1819 PetscFunctionBegin; 1820 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1821 ierr = TSSetUp(ts);CHKERRQ(ierr); 1822 1823 ts->reason = TS_CONVERGED_ITERATING; 1824 1825 ptime_prev = ts->ptime; 1826 ierr = PetscLogEventBegin(TS_Step,ts,0,0,0);CHKERRQ(ierr); 1827 ierr = (*ts->ops->step)(ts);CHKERRQ(ierr); 1828 ierr = PetscLogEventEnd(TS_Step,ts,0,0,0);CHKERRQ(ierr); 1829 ts->time_step_prev = ts->ptime - ptime_prev; 1830 1831 if (ts->reason < 0) { 1832 if (ts->errorifstepfailed) { 1833 if (ts->reason == TS_DIVERGED_NONLINEAR_SOLVE) { 1834 SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s, increase -ts_max_snes_failures or make negative to attempt recovery",TSConvergedReasons[ts->reason]); 1835 } else SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]); 1836 } 1837 } else if (!ts->reason) { 1838 if (ts->steps >= ts->max_steps) 1839 ts->reason = TS_CONVERGED_ITS; 1840 else if (ts->ptime >= ts->max_time) 1841 ts->reason = TS_CONVERGED_TIME; 1842 } 1843 1844 PetscFunctionReturn(0); 1845 } 1846 1847 #undef __FUNCT__ 1848 #define __FUNCT__ "TSEvaluateStep" 1849 /*@ 1850 TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy. 1851 1852 Collective on TS 1853 1854 Input Arguments: 1855 + ts - time stepping context 1856 . order - desired order of accuracy 1857 - done - whether the step was evaluated at this order (pass PETSC_NULL to generate an error if not available) 1858 1859 Output Arguments: 1860 . X - state at the end of the current step 1861 1862 Level: advanced 1863 1864 Notes: 1865 This function cannot be called until all stages have been evaluated. 1866 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. 1867 1868 .seealso: TSStep(), TSAdapt 1869 @*/ 1870 PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec X,PetscBool *done) 1871 { 1872 PetscErrorCode ierr; 1873 1874 PetscFunctionBegin; 1875 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1876 PetscValidType(ts,1); 1877 PetscValidHeaderSpecific(X,VEC_CLASSID,3); 1878 if (!ts->ops->evaluatestep) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSEvaluateStep not implemented for type '%s'",((PetscObject)ts)->type_name); 1879 ierr = (*ts->ops->evaluatestep)(ts,order,X,done);CHKERRQ(ierr); 1880 PetscFunctionReturn(0); 1881 } 1882 1883 #undef __FUNCT__ 1884 #define __FUNCT__ "TSSolve" 1885 /*@ 1886 TSSolve - Steps the requested number of timesteps. 1887 1888 Collective on TS 1889 1890 Input Parameter: 1891 + ts - the TS context obtained from TSCreate() 1892 - x - the solution vector 1893 1894 Output Parameter: 1895 . ftime - time of the state vector x upon completion 1896 1897 Level: beginner 1898 1899 Notes: 1900 The final time returned by this function may be different from the time of the internally 1901 held state accessible by TSGetSolution() and TSGetTime() because the method may have 1902 stepped over the final time. 1903 1904 .keywords: TS, timestep, solve 1905 1906 .seealso: TSCreate(), TSSetSolution(), TSStep() 1907 @*/ 1908 PetscErrorCode TSSolve(TS ts,Vec x,PetscReal *ftime) 1909 { 1910 PetscBool flg; 1911 char filename[PETSC_MAX_PATH_LEN]; 1912 PetscViewer viewer; 1913 PetscErrorCode ierr; 1914 1915 PetscFunctionBegin; 1916 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1917 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 1918 if (ts->exact_final_time) { /* Need ts->vec_sol to be distinct so it is not overwritten when we interpolate at the end */ 1919 if (!ts->vec_sol || x == ts->vec_sol) { 1920 Vec y; 1921 ierr = VecDuplicate(x,&y);CHKERRQ(ierr); 1922 ierr = TSSetSolution(ts,y);CHKERRQ(ierr); 1923 ierr = VecDestroy(&y);CHKERRQ(ierr); /* grant ownership */ 1924 } 1925 ierr = VecCopy(x,ts->vec_sol);CHKERRQ(ierr); 1926 } else { 1927 ierr = TSSetSolution(ts,x);CHKERRQ(ierr); 1928 } 1929 ierr = TSSetUp(ts);CHKERRQ(ierr); 1930 /* reset time step and iteration counters */ 1931 ts->steps = 0; 1932 ts->ksp_its = 0; 1933 ts->snes_its = 0; 1934 ts->num_snes_failures = 0; 1935 ts->reject = 0; 1936 ts->reason = TS_CONVERGED_ITERATING; 1937 1938 if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */ 1939 ierr = (*ts->ops->solve)(ts);CHKERRQ(ierr); 1940 ierr = VecCopy(ts->vec_sol,x);CHKERRQ(ierr); 1941 if (ftime) *ftime = ts->ptime; 1942 } else { 1943 /* steps the requested number of timesteps. */ 1944 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 1945 if (ts->steps >= ts->max_steps) 1946 ts->reason = TS_CONVERGED_ITS; 1947 else if (ts->ptime >= ts->max_time) 1948 ts->reason = TS_CONVERGED_TIME; 1949 while (!ts->reason) { 1950 ierr = TSPreStep(ts);CHKERRQ(ierr); 1951 ierr = TSStep(ts);CHKERRQ(ierr); 1952 ierr = TSPostStep(ts);CHKERRQ(ierr); 1953 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 1954 } 1955 if (ts->exact_final_time && ts->ptime > ts->max_time) { 1956 ierr = TSInterpolate(ts,ts->max_time,x);CHKERRQ(ierr); 1957 if (ftime) *ftime = ts->max_time; 1958 } else { 1959 ierr = VecCopy(ts->vec_sol,x);CHKERRQ(ierr); 1960 if (ftime) *ftime = ts->ptime; 1961 } 1962 } 1963 ierr = PetscOptionsGetString(((PetscObject)ts)->prefix,"-ts_view",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 1964 if (flg && !PetscPreLoadingOn) { 1965 ierr = PetscViewerASCIIOpen(((PetscObject)ts)->comm,filename,&viewer);CHKERRQ(ierr); 1966 ierr = TSView(ts,viewer);CHKERRQ(ierr); 1967 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 1968 } 1969 PetscFunctionReturn(0); 1970 } 1971 1972 #undef __FUNCT__ 1973 #define __FUNCT__ "TSMonitor" 1974 /*@ 1975 TSMonitor - Runs all user-provided monitor routines set using TSMonitorSet() 1976 1977 Collective on TS 1978 1979 Input Parameters: 1980 + ts - time stepping context obtained from TSCreate() 1981 . step - step number that has just completed 1982 . ptime - model time of the state 1983 - x - state at the current model time 1984 1985 Notes: 1986 TSMonitor() is typically used within the time stepping implementations. 1987 Users might call this function when using the TSStep() interface instead of TSSolve(). 1988 1989 Level: advanced 1990 1991 .keywords: TS, timestep 1992 @*/ 1993 PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec x) 1994 { 1995 PetscErrorCode ierr; 1996 PetscInt i,n = ts->numbermonitors; 1997 1998 PetscFunctionBegin; 1999 for (i=0; i<n; i++) { 2000 ierr = (*ts->monitor[i])(ts,step,ptime,x,ts->monitorcontext[i]);CHKERRQ(ierr); 2001 } 2002 PetscFunctionReturn(0); 2003 } 2004 2005 /* ------------------------------------------------------------------------*/ 2006 2007 #undef __FUNCT__ 2008 #define __FUNCT__ "TSMonitorLGCreate" 2009 /*@C 2010 TSMonitorLGCreate - Creates a line graph context for use with 2011 TS to monitor convergence of preconditioned residual norms. 2012 2013 Collective on TS 2014 2015 Input Parameters: 2016 + host - the X display to open, or null for the local machine 2017 . label - the title to put in the title bar 2018 . x, y - the screen coordinates of the upper left coordinate of the window 2019 - m, n - the screen width and height in pixels 2020 2021 Output Parameter: 2022 . draw - the drawing context 2023 2024 Options Database Key: 2025 . -ts_monitor_draw - automatically sets line graph monitor 2026 2027 Notes: 2028 Use TSMonitorLGDestroy() to destroy this line graph, not PetscDrawLGDestroy(). 2029 2030 Level: intermediate 2031 2032 .keywords: TS, monitor, line graph, residual, seealso 2033 2034 .seealso: TSMonitorLGDestroy(), TSMonitorSet() 2035 2036 @*/ 2037 PetscErrorCode TSMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 2038 { 2039 PetscDraw win; 2040 PetscErrorCode ierr; 2041 2042 PetscFunctionBegin; 2043 ierr = PetscDrawCreate(PETSC_COMM_SELF,host,label,x,y,m,n,&win);CHKERRQ(ierr); 2044 ierr = PetscDrawSetType(win,PETSC_DRAW_X);CHKERRQ(ierr); 2045 ierr = PetscDrawLGCreate(win,1,draw);CHKERRQ(ierr); 2046 ierr = PetscDrawLGIndicateDataPoints(*draw);CHKERRQ(ierr); 2047 2048 ierr = PetscLogObjectParent(*draw,win);CHKERRQ(ierr); 2049 PetscFunctionReturn(0); 2050 } 2051 2052 #undef __FUNCT__ 2053 #define __FUNCT__ "TSMonitorLG" 2054 PetscErrorCode TSMonitorLG(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 2055 { 2056 PetscDrawLG lg = (PetscDrawLG) monctx; 2057 PetscReal x,y = ptime; 2058 PetscErrorCode ierr; 2059 2060 PetscFunctionBegin; 2061 if (!monctx) { 2062 MPI_Comm comm; 2063 PetscViewer viewer; 2064 2065 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 2066 viewer = PETSC_VIEWER_DRAW_(comm); 2067 ierr = PetscViewerDrawGetDrawLG(viewer,0,&lg);CHKERRQ(ierr); 2068 } 2069 2070 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 2071 x = (PetscReal)n; 2072 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 2073 if (n < 20 || (n % 5)) { 2074 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 2075 } 2076 PetscFunctionReturn(0); 2077 } 2078 2079 #undef __FUNCT__ 2080 #define __FUNCT__ "TSMonitorLGDestroy" 2081 /*@C 2082 TSMonitorLGDestroy - Destroys a line graph context that was created 2083 with TSMonitorLGCreate(). 2084 2085 Collective on PetscDrawLG 2086 2087 Input Parameter: 2088 . draw - the drawing context 2089 2090 Level: intermediate 2091 2092 .keywords: TS, monitor, line graph, destroy 2093 2094 .seealso: TSMonitorLGCreate(), TSMonitorSet(), TSMonitorLG(); 2095 @*/ 2096 PetscErrorCode TSMonitorLGDestroy(PetscDrawLG *drawlg) 2097 { 2098 PetscDraw draw; 2099 PetscErrorCode ierr; 2100 2101 PetscFunctionBegin; 2102 ierr = PetscDrawLGGetDraw(*drawlg,&draw);CHKERRQ(ierr); 2103 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 2104 ierr = PetscDrawLGDestroy(drawlg);CHKERRQ(ierr); 2105 PetscFunctionReturn(0); 2106 } 2107 2108 #undef __FUNCT__ 2109 #define __FUNCT__ "TSGetTime" 2110 /*@ 2111 TSGetTime - Gets the current time. 2112 2113 Not Collective 2114 2115 Input Parameter: 2116 . ts - the TS context obtained from TSCreate() 2117 2118 Output Parameter: 2119 . t - the current time 2120 2121 Level: beginner 2122 2123 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 2124 2125 .keywords: TS, get, time 2126 @*/ 2127 PetscErrorCode TSGetTime(TS ts,PetscReal* t) 2128 { 2129 PetscFunctionBegin; 2130 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2131 PetscValidDoublePointer(t,2); 2132 *t = ts->ptime; 2133 PetscFunctionReturn(0); 2134 } 2135 2136 #undef __FUNCT__ 2137 #define __FUNCT__ "TSSetTime" 2138 /*@ 2139 TSSetTime - Allows one to reset the time. 2140 2141 Logically Collective on TS 2142 2143 Input Parameters: 2144 + ts - the TS context obtained from TSCreate() 2145 - time - the time 2146 2147 Level: intermediate 2148 2149 .seealso: TSGetTime(), TSSetDuration() 2150 2151 .keywords: TS, set, time 2152 @*/ 2153 PetscErrorCode TSSetTime(TS ts, PetscReal t) 2154 { 2155 PetscFunctionBegin; 2156 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2157 PetscValidLogicalCollectiveReal(ts,t,2); 2158 ts->ptime = t; 2159 PetscFunctionReturn(0); 2160 } 2161 2162 #undef __FUNCT__ 2163 #define __FUNCT__ "TSSetOptionsPrefix" 2164 /*@C 2165 TSSetOptionsPrefix - Sets the prefix used for searching for all 2166 TS options in the database. 2167 2168 Logically Collective on TS 2169 2170 Input Parameter: 2171 + ts - The TS context 2172 - prefix - The prefix to prepend to all option names 2173 2174 Notes: 2175 A hyphen (-) must NOT be given at the beginning of the prefix name. 2176 The first character of all runtime options is AUTOMATICALLY the 2177 hyphen. 2178 2179 Level: advanced 2180 2181 .keywords: TS, set, options, prefix, database 2182 2183 .seealso: TSSetFromOptions() 2184 2185 @*/ 2186 PetscErrorCode TSSetOptionsPrefix(TS ts,const char prefix[]) 2187 { 2188 PetscErrorCode ierr; 2189 SNES snes; 2190 2191 PetscFunctionBegin; 2192 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2193 ierr = PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2194 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2195 ierr = SNESSetOptionsPrefix(snes,prefix);CHKERRQ(ierr); 2196 PetscFunctionReturn(0); 2197 } 2198 2199 2200 #undef __FUNCT__ 2201 #define __FUNCT__ "TSAppendOptionsPrefix" 2202 /*@C 2203 TSAppendOptionsPrefix - Appends to the prefix used for searching for all 2204 TS options in the database. 2205 2206 Logically Collective on TS 2207 2208 Input Parameter: 2209 + ts - The TS context 2210 - prefix - The prefix to prepend to all option names 2211 2212 Notes: 2213 A hyphen (-) must NOT be given at the beginning of the prefix name. 2214 The first character of all runtime options is AUTOMATICALLY the 2215 hyphen. 2216 2217 Level: advanced 2218 2219 .keywords: TS, append, options, prefix, database 2220 2221 .seealso: TSGetOptionsPrefix() 2222 2223 @*/ 2224 PetscErrorCode TSAppendOptionsPrefix(TS ts,const char prefix[]) 2225 { 2226 PetscErrorCode ierr; 2227 SNES snes; 2228 2229 PetscFunctionBegin; 2230 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2231 ierr = PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2232 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2233 ierr = SNESAppendOptionsPrefix(snes,prefix);CHKERRQ(ierr); 2234 PetscFunctionReturn(0); 2235 } 2236 2237 #undef __FUNCT__ 2238 #define __FUNCT__ "TSGetOptionsPrefix" 2239 /*@C 2240 TSGetOptionsPrefix - Sets the prefix used for searching for all 2241 TS options in the database. 2242 2243 Not Collective 2244 2245 Input Parameter: 2246 . ts - The TS context 2247 2248 Output Parameter: 2249 . prefix - A pointer to the prefix string used 2250 2251 Notes: On the fortran side, the user should pass in a string 'prifix' of 2252 sufficient length to hold the prefix. 2253 2254 Level: intermediate 2255 2256 .keywords: TS, get, options, prefix, database 2257 2258 .seealso: TSAppendOptionsPrefix() 2259 @*/ 2260 PetscErrorCode TSGetOptionsPrefix(TS ts,const char *prefix[]) 2261 { 2262 PetscErrorCode ierr; 2263 2264 PetscFunctionBegin; 2265 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2266 PetscValidPointer(prefix,2); 2267 ierr = PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2268 PetscFunctionReturn(0); 2269 } 2270 2271 #undef __FUNCT__ 2272 #define __FUNCT__ "TSGetRHSJacobian" 2273 /*@C 2274 TSGetRHSJacobian - Returns the Jacobian J at the present timestep. 2275 2276 Not Collective, but parallel objects are returned if TS is parallel 2277 2278 Input Parameter: 2279 . ts - The TS context obtained from TSCreate() 2280 2281 Output Parameters: 2282 + J - The Jacobian J of F, where U_t = F(U,t) 2283 . M - The preconditioner matrix, usually the same as J 2284 . func - Function to compute the Jacobian of the RHS 2285 - ctx - User-defined context for Jacobian evaluation routine 2286 2287 Notes: You can pass in PETSC_NULL for any return argument you do not need. 2288 2289 Level: intermediate 2290 2291 .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 2292 2293 .keywords: TS, timestep, get, matrix, Jacobian 2294 @*/ 2295 PetscErrorCode TSGetRHSJacobian(TS ts,Mat *J,Mat *M,TSRHSJacobian *func,void **ctx) 2296 { 2297 PetscErrorCode ierr; 2298 SNES snes; 2299 2300 PetscFunctionBegin; 2301 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2302 ierr = SNESGetJacobian(snes,J,M,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2303 if (func) *func = ts->userops->rhsjacobian; 2304 if (ctx) *ctx = ts->jacP; 2305 PetscFunctionReturn(0); 2306 } 2307 2308 #undef __FUNCT__ 2309 #define __FUNCT__ "TSGetIJacobian" 2310 /*@C 2311 TSGetIJacobian - Returns the implicit Jacobian at the present timestep. 2312 2313 Not Collective, but parallel objects are returned if TS is parallel 2314 2315 Input Parameter: 2316 . ts - The TS context obtained from TSCreate() 2317 2318 Output Parameters: 2319 + A - The Jacobian of F(t,U,U_t) 2320 . B - The preconditioner matrix, often the same as A 2321 . f - The function to compute the matrices 2322 - ctx - User-defined context for Jacobian evaluation routine 2323 2324 Notes: You can pass in PETSC_NULL for any return argument you do not need. 2325 2326 Level: advanced 2327 2328 .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 2329 2330 .keywords: TS, timestep, get, matrix, Jacobian 2331 @*/ 2332 PetscErrorCode TSGetIJacobian(TS ts,Mat *A,Mat *B,TSIJacobian *f,void **ctx) 2333 { 2334 PetscErrorCode ierr; 2335 SNES snes; 2336 2337 PetscFunctionBegin; 2338 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2339 ierr = SNESGetJacobian(snes,A,B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2340 if (f) *f = ts->userops->ijacobian; 2341 if (ctx) *ctx = ts->jacP; 2342 PetscFunctionReturn(0); 2343 } 2344 2345 typedef struct { 2346 PetscViewer viewer; 2347 Vec initialsolution; 2348 PetscBool showinitial; 2349 } TSMonitorSolutionCtx; 2350 2351 #undef __FUNCT__ 2352 #define __FUNCT__ "TSMonitorSolution" 2353 /*@C 2354 TSMonitorSolution - Monitors progress of the TS solvers by calling 2355 VecView() for the solution at each timestep 2356 2357 Collective on TS 2358 2359 Input Parameters: 2360 + ts - the TS context 2361 . step - current time-step 2362 . ptime - current time 2363 - dummy - either a viewer or PETSC_NULL 2364 2365 Level: intermediate 2366 2367 .keywords: TS, vector, monitor, view 2368 2369 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 2370 @*/ 2371 PetscErrorCode TSMonitorSolution(TS ts,PetscInt step,PetscReal ptime,Vec x,void *dummy) 2372 { 2373 PetscErrorCode ierr; 2374 TSMonitorSolutionCtx *ictx = (TSMonitorSolutionCtx*)dummy; 2375 2376 PetscFunctionBegin; 2377 if (!step && ictx->showinitial) { 2378 if (!ictx->initialsolution) { 2379 ierr = VecDuplicate(x,&ictx->initialsolution);CHKERRQ(ierr); 2380 } 2381 ierr = VecCopy(x,ictx->initialsolution);CHKERRQ(ierr); 2382 } 2383 if (ictx->showinitial) { 2384 PetscReal pause; 2385 ierr = PetscViewerDrawGetPause(ictx->viewer,&pause);CHKERRQ(ierr); 2386 ierr = PetscViewerDrawSetPause(ictx->viewer,0.0);CHKERRQ(ierr); 2387 ierr = VecView(ictx->initialsolution,ictx->viewer);CHKERRQ(ierr); 2388 ierr = PetscViewerDrawSetPause(ictx->viewer,pause);CHKERRQ(ierr); 2389 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);CHKERRQ(ierr); 2390 } 2391 ierr = VecView(x,ictx->viewer);CHKERRQ(ierr); 2392 if (ictx->showinitial) { 2393 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);CHKERRQ(ierr); 2394 } 2395 PetscFunctionReturn(0); 2396 } 2397 2398 2399 #undef __FUNCT__ 2400 #define __FUNCT__ "TSMonitorSolutionDestroy" 2401 /*@C 2402 TSMonitorSolutionDestroy - Destroys the monitor context for TSMonitorSolution 2403 2404 Collective on TS 2405 2406 Input Parameters: 2407 . ctx - the monitor context 2408 2409 Level: intermediate 2410 2411 .keywords: TS, vector, monitor, view 2412 2413 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorSolution() 2414 @*/ 2415 PetscErrorCode TSMonitorSolutionDestroy(void **ctx) 2416 { 2417 PetscErrorCode ierr; 2418 TSMonitorSolutionCtx *ictx = *(TSMonitorSolutionCtx**)ctx; 2419 2420 PetscFunctionBegin; 2421 ierr = PetscViewerDestroy(&ictx->viewer);CHKERRQ(ierr); 2422 ierr = VecDestroy(&ictx->initialsolution);CHKERRQ(ierr); 2423 ierr = PetscFree(ictx);CHKERRQ(ierr); 2424 PetscFunctionReturn(0); 2425 } 2426 2427 #undef __FUNCT__ 2428 #define __FUNCT__ "TSMonitorSolutionCreate" 2429 /*@C 2430 TSMonitorSolutionCreate - Creates the monitor context for TSMonitorSolution 2431 2432 Collective on TS 2433 2434 Input Parameter: 2435 . ts - time-step context 2436 2437 Output Patameter: 2438 . ctx - the monitor context 2439 2440 Level: intermediate 2441 2442 .keywords: TS, vector, monitor, view 2443 2444 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorSolution() 2445 @*/ 2446 PetscErrorCode TSMonitorSolutionCreate(TS ts,PetscViewer viewer,void **ctx) 2447 { 2448 PetscErrorCode ierr; 2449 TSMonitorSolutionCtx *ictx; 2450 2451 PetscFunctionBegin; 2452 ierr = PetscNew(TSMonitorSolutionCtx,&ictx);CHKERRQ(ierr); 2453 *ctx = (void*)ictx; 2454 if (!viewer) { 2455 viewer = PETSC_VIEWER_DRAW_(((PetscObject)ts)->comm); 2456 } 2457 ierr = PetscObjectReference((PetscObject)viewer);CHKERRQ(ierr); 2458 ictx->viewer = viewer; 2459 ictx->showinitial = PETSC_FALSE; 2460 ierr = PetscOptionsGetBool(((PetscObject)ts)->prefix,"-ts_monitor_solution_initial",&ictx->showinitial,PETSC_NULL);CHKERRQ(ierr); 2461 PetscFunctionReturn(0); 2462 } 2463 2464 #undef __FUNCT__ 2465 #define __FUNCT__ "TSSetDM" 2466 /*@ 2467 TSSetDM - Sets the DM that may be used by some preconditioners 2468 2469 Logically Collective on TS and DM 2470 2471 Input Parameters: 2472 + ts - the preconditioner context 2473 - dm - the dm 2474 2475 Level: intermediate 2476 2477 2478 .seealso: TSGetDM(), SNESSetDM(), SNESGetDM() 2479 @*/ 2480 PetscErrorCode TSSetDM(TS ts,DM dm) 2481 { 2482 PetscErrorCode ierr; 2483 SNES snes; 2484 2485 PetscFunctionBegin; 2486 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2487 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 2488 ierr = DMDestroy(&ts->dm);CHKERRQ(ierr); 2489 ts->dm = dm; 2490 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2491 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 2492 PetscFunctionReturn(0); 2493 } 2494 2495 #undef __FUNCT__ 2496 #define __FUNCT__ "TSGetDM" 2497 /*@ 2498 TSGetDM - Gets the DM that may be used by some preconditioners 2499 2500 Not Collective 2501 2502 Input Parameter: 2503 . ts - the preconditioner context 2504 2505 Output Parameter: 2506 . dm - the dm 2507 2508 Level: intermediate 2509 2510 2511 .seealso: TSSetDM(), SNESSetDM(), SNESGetDM() 2512 @*/ 2513 PetscErrorCode TSGetDM(TS ts,DM *dm) 2514 { 2515 PetscErrorCode ierr; 2516 2517 PetscFunctionBegin; 2518 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2519 if (!ts->dm) { 2520 ierr = DMShellCreate(((PetscObject)ts)->comm,&ts->dm);CHKERRQ(ierr); 2521 if (ts->snes) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 2522 } 2523 *dm = ts->dm; 2524 PetscFunctionReturn(0); 2525 } 2526 2527 #undef __FUNCT__ 2528 #define __FUNCT__ "SNESTSFormFunction" 2529 /*@ 2530 SNESTSFormFunction - Function to evaluate nonlinear residual 2531 2532 Logically Collective on SNES 2533 2534 Input Parameter: 2535 + snes - nonlinear solver 2536 . X - the current state at which to evaluate the residual 2537 - ctx - user context, must be a TS 2538 2539 Output Parameter: 2540 . F - the nonlinear residual 2541 2542 Notes: 2543 This function is not normally called by users and is automatically registered with the SNES used by TS. 2544 It is most frequently passed to MatFDColoringSetFunction(). 2545 2546 Level: advanced 2547 2548 .seealso: SNESSetFunction(), MatFDColoringSetFunction() 2549 @*/ 2550 PetscErrorCode SNESTSFormFunction(SNES snes,Vec X,Vec F,void *ctx) 2551 { 2552 TS ts = (TS)ctx; 2553 PetscErrorCode ierr; 2554 2555 PetscFunctionBegin; 2556 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2557 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2558 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 2559 PetscValidHeaderSpecific(ts,TS_CLASSID,4); 2560 ierr = (ts->ops->snesfunction)(snes,X,F,ts);CHKERRQ(ierr); 2561 PetscFunctionReturn(0); 2562 } 2563 2564 #undef __FUNCT__ 2565 #define __FUNCT__ "SNESTSFormJacobian" 2566 /*@ 2567 SNESTSFormJacobian - Function to evaluate the Jacobian 2568 2569 Collective on SNES 2570 2571 Input Parameter: 2572 + snes - nonlinear solver 2573 . X - the current state at which to evaluate the residual 2574 - ctx - user context, must be a TS 2575 2576 Output Parameter: 2577 + A - the Jacobian 2578 . B - the preconditioning matrix (may be the same as A) 2579 - flag - indicates any structure change in the matrix 2580 2581 Notes: 2582 This function is not normally called by users and is automatically registered with the SNES used by TS. 2583 2584 Level: developer 2585 2586 .seealso: SNESSetJacobian() 2587 @*/ 2588 PetscErrorCode SNESTSFormJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flag,void *ctx) 2589 { 2590 TS ts = (TS)ctx; 2591 PetscErrorCode ierr; 2592 2593 PetscFunctionBegin; 2594 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2595 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2596 PetscValidPointer(A,3); 2597 PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 2598 PetscValidPointer(B,4); 2599 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); 2600 PetscValidPointer(flag,5); 2601 PetscValidHeaderSpecific(ts,TS_CLASSID,6); 2602 ierr = (ts->ops->snesjacobian)(snes,X,A,B,flag,ts);CHKERRQ(ierr); 2603 PetscFunctionReturn(0); 2604 } 2605 2606 #undef __FUNCT__ 2607 #define __FUNCT__ "TSComputeRHSFunctionLinear" 2608 /*@C 2609 TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems only 2610 2611 Collective on TS 2612 2613 Input Arguments: 2614 + ts - time stepping context 2615 . t - time at which to evaluate 2616 . X - state at which to evaluate 2617 - ctx - context 2618 2619 Output Arguments: 2620 . F - right hand side 2621 2622 Level: intermediate 2623 2624 Notes: 2625 This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems. 2626 The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian(). 2627 2628 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant() 2629 @*/ 2630 PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec X,Vec F,void *ctx) 2631 { 2632 PetscErrorCode ierr; 2633 Mat Arhs,Brhs; 2634 MatStructure flg2; 2635 2636 PetscFunctionBegin; 2637 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 2638 ierr = TSComputeRHSJacobian(ts,t,X,&Arhs,&Brhs,&flg2);CHKERRQ(ierr); 2639 ierr = MatMult(Arhs,X,F);CHKERRQ(ierr); 2640 PetscFunctionReturn(0); 2641 } 2642 2643 #undef __FUNCT__ 2644 #define __FUNCT__ "TSComputeRHSJacobianConstant" 2645 /*@C 2646 TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent. 2647 2648 Collective on TS 2649 2650 Input Arguments: 2651 + ts - time stepping context 2652 . t - time at which to evaluate 2653 . X - state at which to evaluate 2654 - ctx - context 2655 2656 Output Arguments: 2657 + A - pointer to operator 2658 . B - pointer to preconditioning matrix 2659 - flg - matrix structure flag 2660 2661 Level: intermediate 2662 2663 Notes: 2664 This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems. 2665 2666 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear() 2667 @*/ 2668 PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec X,Mat *A,Mat *B,MatStructure *flg,void *ctx) 2669 { 2670 2671 PetscFunctionBegin; 2672 *flg = SAME_PRECONDITIONER; 2673 PetscFunctionReturn(0); 2674 } 2675 2676 #undef __FUNCT__ 2677 #define __FUNCT__ "TSComputeIFunctionLinear" 2678 /*@C 2679 TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only 2680 2681 Collective on TS 2682 2683 Input Arguments: 2684 + ts - time stepping context 2685 . t - time at which to evaluate 2686 . X - state at which to evaluate 2687 . Xdot - time derivative of state vector 2688 - ctx - context 2689 2690 Output Arguments: 2691 . F - left hand side 2692 2693 Level: intermediate 2694 2695 Notes: 2696 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 2697 user is required to write their own TSComputeIFunction. 2698 This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems. 2699 The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian(). 2700 2701 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant() 2702 @*/ 2703 PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx) 2704 { 2705 PetscErrorCode ierr; 2706 Mat A,B; 2707 MatStructure flg2; 2708 2709 PetscFunctionBegin; 2710 ierr = TSGetIJacobian(ts,&A,&B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2711 ierr = TSComputeIJacobian(ts,t,X,Xdot,1.0,&A,&B,&flg2,PETSC_TRUE);CHKERRQ(ierr); 2712 ierr = MatMult(A,Xdot,F);CHKERRQ(ierr); 2713 PetscFunctionReturn(0); 2714 } 2715 2716 #undef __FUNCT__ 2717 #define __FUNCT__ "TSComputeIJacobianConstant" 2718 /*@C 2719 TSComputeIJacobianConstant - Reuses a Jacobian that is time-independent. 2720 2721 Collective on TS 2722 2723 Input Arguments: 2724 + ts - time stepping context 2725 . t - time at which to evaluate 2726 . X - state at which to evaluate 2727 . Xdot - time derivative of state vector 2728 . shift - shift to apply 2729 - ctx - context 2730 2731 Output Arguments: 2732 + A - pointer to operator 2733 . B - pointer to preconditioning matrix 2734 - flg - matrix structure flag 2735 2736 Level: intermediate 2737 2738 Notes: 2739 This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems. 2740 2741 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear() 2742 @*/ 2743 PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,void *ctx) 2744 { 2745 2746 PetscFunctionBegin; 2747 *flg = SAME_PRECONDITIONER; 2748 PetscFunctionReturn(0); 2749 } 2750 2751 2752 #undef __FUNCT__ 2753 #define __FUNCT__ "TSGetConvergedReason" 2754 /*@ 2755 TSGetConvergedReason - Gets the reason the TS iteration was stopped. 2756 2757 Not Collective 2758 2759 Input Parameter: 2760 . ts - the TS context 2761 2762 Output Parameter: 2763 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 2764 manual pages for the individual convergence tests for complete lists 2765 2766 Level: intermediate 2767 2768 Notes: 2769 Can only be called after the call to TSSolve() is complete. 2770 2771 .keywords: TS, nonlinear, set, convergence, test 2772 2773 .seealso: TSSetConvergenceTest(), TSConvergedReason 2774 @*/ 2775 PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason) 2776 { 2777 PetscFunctionBegin; 2778 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2779 PetscValidPointer(reason,2); 2780 *reason = ts->reason; 2781 PetscFunctionReturn(0); 2782 } 2783 2784 #undef __FUNCT__ 2785 #define __FUNCT__ "TSGetSNESIterations" 2786 /*@ 2787 TSGetSNESIterations - Gets the total number of nonlinear iterations 2788 used by the time integrator. 2789 2790 Not Collective 2791 2792 Input Parameter: 2793 . ts - TS context 2794 2795 Output Parameter: 2796 . nits - number of nonlinear iterations 2797 2798 Notes: 2799 This counter is reset to zero for each successive call to TSSolve(). 2800 2801 Level: intermediate 2802 2803 .keywords: TS, get, number, nonlinear, iterations 2804 2805 .seealso: TSGetKSPIterations() 2806 @*/ 2807 PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits) 2808 { 2809 PetscFunctionBegin; 2810 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2811 PetscValidIntPointer(nits,2); 2812 *nits = ts->snes_its; 2813 PetscFunctionReturn(0); 2814 } 2815 2816 #undef __FUNCT__ 2817 #define __FUNCT__ "TSGetKSPIterations" 2818 /*@ 2819 TSGetKSPIterations - Gets the total number of linear iterations 2820 used by the time integrator. 2821 2822 Not Collective 2823 2824 Input Parameter: 2825 . ts - TS context 2826 2827 Output Parameter: 2828 . lits - number of linear iterations 2829 2830 Notes: 2831 This counter is reset to zero for each successive call to TSSolve(). 2832 2833 Level: intermediate 2834 2835 .keywords: TS, get, number, linear, iterations 2836 2837 .seealso: TSGetSNESIterations(), SNESGetKSPIterations() 2838 @*/ 2839 PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits) 2840 { 2841 PetscFunctionBegin; 2842 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2843 PetscValidIntPointer(lits,2); 2844 *lits = ts->ksp_its; 2845 PetscFunctionReturn(0); 2846 } 2847 2848 #undef __FUNCT__ 2849 #define __FUNCT__ "TSGetStepRejections" 2850 /*@ 2851 TSGetStepRejections - Gets the total number of rejected steps. 2852 2853 Not Collective 2854 2855 Input Parameter: 2856 . ts - TS context 2857 2858 Output Parameter: 2859 . rejects - number of steps rejected 2860 2861 Notes: 2862 This counter is reset to zero for each successive call to TSSolve(). 2863 2864 Level: intermediate 2865 2866 .keywords: TS, get, number 2867 2868 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails() 2869 @*/ 2870 PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects) 2871 { 2872 PetscFunctionBegin; 2873 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2874 PetscValidIntPointer(rejects,2); 2875 *rejects = ts->reject; 2876 PetscFunctionReturn(0); 2877 } 2878 2879 #undef __FUNCT__ 2880 #define __FUNCT__ "TSGetSNESFailures" 2881 /*@ 2882 TSGetSNESFailures - Gets the total number of failed SNES solves 2883 2884 Not Collective 2885 2886 Input Parameter: 2887 . ts - TS context 2888 2889 Output Parameter: 2890 . fails - number of failed nonlinear solves 2891 2892 Notes: 2893 This counter is reset to zero for each successive call to TSSolve(). 2894 2895 Level: intermediate 2896 2897 .keywords: TS, get, number 2898 2899 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures() 2900 @*/ 2901 PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails) 2902 { 2903 PetscFunctionBegin; 2904 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2905 PetscValidIntPointer(fails,2); 2906 *fails = ts->num_snes_failures; 2907 PetscFunctionReturn(0); 2908 } 2909 2910 #undef __FUNCT__ 2911 #define __FUNCT__ "TSSetMaxStepRejections" 2912 /*@ 2913 TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails 2914 2915 Not Collective 2916 2917 Input Parameter: 2918 + ts - TS context 2919 - rejects - maximum number of rejected steps, pass -1 for unlimited 2920 2921 Notes: 2922 The counter is reset to zero for each step 2923 2924 Options Database Key: 2925 . -ts_max_reject - Maximum number of step rejections before a step fails 2926 2927 Level: intermediate 2928 2929 .keywords: TS, set, maximum, number 2930 2931 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 2932 @*/ 2933 PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects) 2934 { 2935 PetscFunctionBegin; 2936 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2937 ts->max_reject = rejects; 2938 PetscFunctionReturn(0); 2939 } 2940 2941 #undef __FUNCT__ 2942 #define __FUNCT__ "TSSetMaxSNESFailures" 2943 /*@ 2944 TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves 2945 2946 Not Collective 2947 2948 Input Parameter: 2949 + ts - TS context 2950 - fails - maximum number of failed nonlinear solves, pass -1 for unlimited 2951 2952 Notes: 2953 The counter is reset to zero for each successive call to TSSolve(). 2954 2955 Options Database Key: 2956 . -ts_max_snes_failures - Maximum number of nonlinear solve failures 2957 2958 Level: intermediate 2959 2960 .keywords: TS, set, maximum, number 2961 2962 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason() 2963 @*/ 2964 PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails) 2965 { 2966 PetscFunctionBegin; 2967 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2968 ts->max_snes_failures = fails; 2969 PetscFunctionReturn(0); 2970 } 2971 2972 #undef __FUNCT__ 2973 #define __FUNCT__ "TSSetErrorIfStepFails()" 2974 /*@ 2975 TSSetErrorIfStepFails - Error if no step succeeds 2976 2977 Not Collective 2978 2979 Input Parameter: 2980 + ts - TS context 2981 - err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure 2982 2983 Options Database Key: 2984 . -ts_error_if_step_fails - Error if no step succeeds 2985 2986 Level: intermediate 2987 2988 .keywords: TS, set, error 2989 2990 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 2991 @*/ 2992 PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err) 2993 { 2994 PetscFunctionBegin; 2995 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2996 ts->errorifstepfailed = err; 2997 PetscFunctionReturn(0); 2998 } 2999 3000 #undef __FUNCT__ 3001 #define __FUNCT__ "TSMonitorSolutionBinary" 3002 /*@C 3003 TSMonitorSolutionBinary - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file 3004 3005 Collective on TS 3006 3007 Input Parameters: 3008 + ts - the TS context 3009 . step - current time-step 3010 . ptime - current time 3011 . x - current state 3012 - viewer - binary viewer 3013 3014 Level: intermediate 3015 3016 .keywords: TS, vector, monitor, view 3017 3018 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3019 @*/ 3020 PetscErrorCode TSMonitorSolutionBinary(TS ts,PetscInt step,PetscReal ptime,Vec x,void *viewer) 3021 { 3022 PetscErrorCode ierr; 3023 PetscViewer v = (PetscViewer)viewer; 3024 3025 PetscFunctionBegin; 3026 ierr = VecView(x,v);CHKERRQ(ierr); 3027 PetscFunctionReturn(0); 3028 } 3029 3030 #undef __FUNCT__ 3031 #define __FUNCT__ "TSMonitorSolutionVTK" 3032 /*@C 3033 TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep. 3034 3035 Collective on TS 3036 3037 Input Parameters: 3038 + ts - the TS context 3039 . step - current time-step 3040 . ptime - current time 3041 . x - current state 3042 - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 3043 3044 Level: intermediate 3045 3046 Notes: 3047 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. 3048 These are named according to the file name template. 3049 3050 This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy(). 3051 3052 .keywords: TS, vector, monitor, view 3053 3054 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3055 @*/ 3056 PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec x,void *filenametemplate) 3057 { 3058 PetscErrorCode ierr; 3059 char filename[PETSC_MAX_PATH_LEN]; 3060 PetscViewer viewer; 3061 3062 PetscFunctionBegin; 3063 ierr = PetscSNPrintf(filename,sizeof filename,(const char*)filenametemplate,step);CHKERRQ(ierr); 3064 ierr = PetscViewerVTKOpen(((PetscObject)ts)->comm,filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 3065 ierr = VecView(x,viewer);CHKERRQ(ierr); 3066 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3067 PetscFunctionReturn(0); 3068 } 3069 3070 #undef __FUNCT__ 3071 #define __FUNCT__ "TSMonitorSolutionVTKDestroy" 3072 /*@C 3073 TSMonitorSolutionVTKDestroy - Destroy context for monitoring 3074 3075 Collective on TS 3076 3077 Input Parameters: 3078 . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 3079 3080 Level: intermediate 3081 3082 Note: 3083 This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK(). 3084 3085 .keywords: TS, vector, monitor, view 3086 3087 .seealso: TSMonitorSet(), TSMonitorSolutionVTK() 3088 @*/ 3089 PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate) 3090 { 3091 PetscErrorCode ierr; 3092 3093 PetscFunctionBegin; 3094 ierr = PetscFree(*(char**)filenametemplate);CHKERRQ(ierr); 3095 PetscFunctionReturn(0); 3096 } 3097 3098 #undef __FUNCT__ 3099 #define __FUNCT__ "TSGetAdapt" 3100 /*@ 3101 TSGetAdapt - Get the adaptive controller context for the current method 3102 3103 Collective on TS if controller has not been created yet 3104 3105 Input Arguments: 3106 . ts - time stepping context 3107 3108 Output Arguments: 3109 . adapt - adaptive controller 3110 3111 Level: intermediate 3112 3113 .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose() 3114 @*/ 3115 PetscErrorCode TSGetAdapt(TS ts,TSAdapt *adapt) 3116 { 3117 PetscErrorCode ierr; 3118 3119 PetscFunctionBegin; 3120 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3121 PetscValidPointer(adapt,2); 3122 if (!ts->adapt) { 3123 ierr = TSAdaptCreate(((PetscObject)ts)->comm,&ts->adapt);CHKERRQ(ierr); 3124 ierr = PetscLogObjectParent(ts,ts->adapt);CHKERRQ(ierr); 3125 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);CHKERRQ(ierr); 3126 } 3127 *adapt = ts->adapt; 3128 PetscFunctionReturn(0); 3129 } 3130 3131 #undef __FUNCT__ 3132 #define __FUNCT__ "TSSetTolerances" 3133 /*@ 3134 TSSetTolerances - Set tolerances for local truncation error when using adaptive controller 3135 3136 Logically Collective 3137 3138 Input Arguments: 3139 + ts - time integration context 3140 . atol - scalar absolute tolerances, PETSC_DECIDE to leave current value 3141 . vatol - vector of absolute tolerances or PETSC_NULL, used in preference to atol if present 3142 . rtol - scalar relative tolerances, PETSC_DECIDE to leave current value 3143 - vrtol - vector of relative tolerances or PETSC_NULL, used in preference to atol if present 3144 3145 Level: beginner 3146 3147 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSGetTolerances() 3148 @*/ 3149 PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol) 3150 { 3151 PetscErrorCode ierr; 3152 3153 PetscFunctionBegin; 3154 if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol; 3155 if (vatol) { 3156 ierr = PetscObjectReference((PetscObject)vatol);CHKERRQ(ierr); 3157 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 3158 ts->vatol = vatol; 3159 } 3160 if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol; 3161 if (vrtol) { 3162 ierr = PetscObjectReference((PetscObject)vrtol);CHKERRQ(ierr); 3163 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 3164 ts->vrtol = vrtol; 3165 } 3166 PetscFunctionReturn(0); 3167 } 3168 3169 #undef __FUNCT__ 3170 #define __FUNCT__ "TSGetTolerances" 3171 /*@ 3172 TSGetTolerances - Get tolerances for local truncation error when using adaptive controller 3173 3174 Logically Collective 3175 3176 Input Arguments: 3177 . ts - time integration context 3178 3179 Output Arguments: 3180 + atol - scalar absolute tolerances, PETSC_NULL to ignore 3181 . vatol - vector of absolute tolerances, PETSC_NULL to ignore 3182 . rtol - scalar relative tolerances, PETSC_NULL to ignore 3183 - vrtol - vector of relative tolerances, PETSC_NULL to ignore 3184 3185 Level: beginner 3186 3187 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSSetTolerances() 3188 @*/ 3189 PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol) 3190 { 3191 3192 PetscFunctionBegin; 3193 if (atol) *atol = ts->atol; 3194 if (vatol) *vatol = ts->vatol; 3195 if (rtol) *rtol = ts->rtol; 3196 if (vrtol) *vrtol = ts->vrtol; 3197 PetscFunctionReturn(0); 3198 } 3199 3200 #undef __FUNCT__ 3201 #define __FUNCT__ "TSErrorNormWRMS" 3202 /*@ 3203 TSErrorNormWRMS - compute a weighted norm of the difference between a vector and the current state 3204 3205 Collective on TS 3206 3207 Input Arguments: 3208 + ts - time stepping context 3209 - Y - state vector to be compared to ts->vec_sol 3210 3211 Output Arguments: 3212 . norm - weighted norm, a value of 1.0 is considered small 3213 3214 Level: developer 3215 3216 .seealso: TSSetTolerances() 3217 @*/ 3218 PetscErrorCode TSErrorNormWRMS(TS ts,Vec Y,PetscReal *norm) 3219 { 3220 PetscErrorCode ierr; 3221 PetscInt i,n,N; 3222 const PetscScalar *x,*y; 3223 Vec X; 3224 PetscReal sum,gsum; 3225 3226 PetscFunctionBegin; 3227 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3228 PetscValidHeaderSpecific(Y,VEC_CLASSID,2); 3229 PetscValidPointer(norm,3); 3230 X = ts->vec_sol; 3231 PetscCheckSameTypeAndComm(X,1,Y,2); 3232 if (X == Y) SETERRQ(((PetscObject)X)->comm,PETSC_ERR_ARG_IDN,"Y cannot be the TS solution vector"); 3233 3234 ierr = VecGetSize(X,&N);CHKERRQ(ierr); 3235 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr); 3236 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 3237 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 3238 sum = 0.; 3239 if (ts->vatol && ts->vrtol) { 3240 const PetscScalar *atol,*rtol; 3241 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3242 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3243 for (i=0; i<n; i++) { 3244 PetscReal tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(x[i]),PetscAbsScalar(y[i])); 3245 sum += PetscSqr(PetscAbsScalar(y[i] - x[i]) / tol); 3246 } 3247 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3248 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3249 } else if (ts->vatol) { /* vector atol, scalar rtol */ 3250 const PetscScalar *atol; 3251 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3252 for (i=0; i<n; i++) { 3253 PetscReal tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(x[i]),PetscAbsScalar(y[i])); 3254 sum += PetscSqr(PetscAbsScalar(y[i] - x[i]) / tol); 3255 } 3256 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3257 } else if (ts->vrtol) { /* scalar atol, vector rtol */ 3258 const PetscScalar *rtol; 3259 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3260 for (i=0; i<n; i++) { 3261 PetscReal tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(x[i]),PetscAbsScalar(y[i])); 3262 sum += PetscSqr(PetscAbsScalar(y[i] - x[i]) / tol); 3263 } 3264 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3265 } else { /* scalar atol, scalar rtol */ 3266 for (i=0; i<n; i++) { 3267 PetscReal tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(x[i]),PetscAbsScalar(y[i])); 3268 sum += PetscSqr(PetscAbsScalar(y[i] - x[i]) / tol); 3269 } 3270 } 3271 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 3272 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 3273 3274 ierr = MPI_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,((PetscObject)ts)->comm);CHKERRQ(ierr); 3275 *norm = PetscSqrtReal(gsum / N); 3276 if (PetscIsInfOrNanScalar(*norm)) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 3277 PetscFunctionReturn(0); 3278 } 3279 3280 #undef __FUNCT__ 3281 #define __FUNCT__ "TSSetCFLTimeLocal" 3282 /*@ 3283 TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler 3284 3285 Logically Collective on TS 3286 3287 Input Arguments: 3288 + ts - time stepping context 3289 - cfltime - maximum stable time step if using forward Euler (value can be different on each process) 3290 3291 Note: 3292 After calling this function, the global CFL time can be obtained by calling TSGetCFLTime() 3293 3294 Level: intermediate 3295 3296 .seealso: TSGetCFLTime(), TSADAPTCFL 3297 @*/ 3298 PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime) 3299 { 3300 3301 PetscFunctionBegin; 3302 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3303 ts->cfltime_local = cfltime; 3304 ts->cfltime = -1.; 3305 PetscFunctionReturn(0); 3306 } 3307 3308 #undef __FUNCT__ 3309 #define __FUNCT__ "TSGetCFLTime" 3310 /*@ 3311 TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler 3312 3313 Collective on TS 3314 3315 Input Arguments: 3316 . ts - time stepping context 3317 3318 Output Arguments: 3319 . cfltime - maximum stable time step for forward Euler 3320 3321 Level: advanced 3322 3323 .seealso: TSSetCFLTimeLocal() 3324 @*/ 3325 PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime) 3326 { 3327 PetscErrorCode ierr; 3328 3329 PetscFunctionBegin; 3330 if (ts->cfltime < 0) { 3331 ierr = MPI_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,((PetscObject)ts)->comm);CHKERRQ(ierr); 3332 } 3333 *cfltime = ts->cfltime; 3334 PetscFunctionReturn(0); 3335 } 3336 3337 #undef __FUNCT__ 3338 #define __FUNCT__ "TSVISetVariableBounds" 3339 /*@ 3340 TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu 3341 3342 Input Parameters: 3343 . ts - the TS context. 3344 . xl - lower bound. 3345 . xu - upper bound. 3346 3347 Notes: 3348 If this routine is not called then the lower and upper bounds are set to 3349 SNES_VI_NINF and SNES_VI_INF respectively during SNESSetUp(). 3350 3351 Level: advanced 3352 3353 @*/ 3354 PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu) 3355 { 3356 PetscErrorCode ierr; 3357 SNES snes; 3358 3359 PetscFunctionBegin; 3360 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3361 ierr = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(ierr); 3362 PetscFunctionReturn(0); 3363 } 3364 3365 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3366 #include <mex.h> 3367 3368 typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext; 3369 3370 #undef __FUNCT__ 3371 #define __FUNCT__ "TSComputeFunction_Matlab" 3372 /* 3373 TSComputeFunction_Matlab - Calls the function that has been set with 3374 TSSetFunctionMatlab(). 3375 3376 Collective on TS 3377 3378 Input Parameters: 3379 + snes - the TS context 3380 - x - input vector 3381 3382 Output Parameter: 3383 . y - function vector, as set by TSSetFunction() 3384 3385 Notes: 3386 TSComputeFunction() is typically used within nonlinear solvers 3387 implementations, so most users would not generally call this routine 3388 themselves. 3389 3390 Level: developer 3391 3392 .keywords: TS, nonlinear, compute, function 3393 3394 .seealso: TSSetFunction(), TSGetFunction() 3395 */ 3396 PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec x,Vec xdot,Vec y, void *ctx) 3397 { 3398 PetscErrorCode ierr; 3399 TSMatlabContext *sctx = (TSMatlabContext *)ctx; 3400 int nlhs = 1,nrhs = 7; 3401 mxArray *plhs[1],*prhs[7]; 3402 long long int lx = 0,lxdot = 0,ly = 0,ls = 0; 3403 3404 PetscFunctionBegin; 3405 PetscValidHeaderSpecific(snes,TS_CLASSID,1); 3406 PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3407 PetscValidHeaderSpecific(xdot,VEC_CLASSID,4); 3408 PetscValidHeaderSpecific(y,VEC_CLASSID,5); 3409 PetscCheckSameComm(snes,1,x,3); 3410 PetscCheckSameComm(snes,1,y,5); 3411 3412 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 3413 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 3414 ierr = PetscMemcpy(&lxdot,&xdot,sizeof(xdot));CHKERRQ(ierr); 3415 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 3416 prhs[0] = mxCreateDoubleScalar((double)ls); 3417 prhs[1] = mxCreateDoubleScalar(time); 3418 prhs[2] = mxCreateDoubleScalar((double)lx); 3419 prhs[3] = mxCreateDoubleScalar((double)lxdot); 3420 prhs[4] = mxCreateDoubleScalar((double)ly); 3421 prhs[5] = mxCreateString(sctx->funcname); 3422 prhs[6] = sctx->ctx; 3423 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");CHKERRQ(ierr); 3424 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 3425 mxDestroyArray(prhs[0]); 3426 mxDestroyArray(prhs[1]); 3427 mxDestroyArray(prhs[2]); 3428 mxDestroyArray(prhs[3]); 3429 mxDestroyArray(prhs[4]); 3430 mxDestroyArray(prhs[5]); 3431 mxDestroyArray(plhs[0]); 3432 PetscFunctionReturn(0); 3433 } 3434 3435 3436 #undef __FUNCT__ 3437 #define __FUNCT__ "TSSetFunctionMatlab" 3438 /* 3439 TSSetFunctionMatlab - Sets the function evaluation routine and function 3440 vector for use by the TS routines in solving ODEs 3441 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 3442 3443 Logically Collective on TS 3444 3445 Input Parameters: 3446 + ts - the TS context 3447 - func - function evaluation routine 3448 3449 Calling sequence of func: 3450 $ func (TS ts,PetscReal time,Vec x,Vec xdot,Vec f,void *ctx); 3451 3452 Level: beginner 3453 3454 .keywords: TS, nonlinear, set, function 3455 3456 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 3457 */ 3458 PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx) 3459 { 3460 PetscErrorCode ierr; 3461 TSMatlabContext *sctx; 3462 3463 PetscFunctionBegin; 3464 /* currently sctx is memory bleed */ 3465 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 3466 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 3467 /* 3468 This should work, but it doesn't 3469 sctx->ctx = ctx; 3470 mexMakeArrayPersistent(sctx->ctx); 3471 */ 3472 sctx->ctx = mxDuplicateArray(ctx); 3473 ierr = TSSetIFunction(ts,PETSC_NULL,TSComputeFunction_Matlab,sctx);CHKERRQ(ierr); 3474 PetscFunctionReturn(0); 3475 } 3476 3477 #undef __FUNCT__ 3478 #define __FUNCT__ "TSComputeJacobian_Matlab" 3479 /* 3480 TSComputeJacobian_Matlab - Calls the function that has been set with 3481 TSSetJacobianMatlab(). 3482 3483 Collective on TS 3484 3485 Input Parameters: 3486 + ts - the TS context 3487 . x - input vector 3488 . A, B - the matrices 3489 - ctx - user context 3490 3491 Output Parameter: 3492 . flag - structure of the matrix 3493 3494 Level: developer 3495 3496 .keywords: TS, nonlinear, compute, function 3497 3498 .seealso: TSSetFunction(), TSGetFunction() 3499 @*/ 3500 PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec x,Vec xdot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag, void *ctx) 3501 { 3502 PetscErrorCode ierr; 3503 TSMatlabContext *sctx = (TSMatlabContext *)ctx; 3504 int nlhs = 2,nrhs = 9; 3505 mxArray *plhs[2],*prhs[9]; 3506 long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0; 3507 3508 PetscFunctionBegin; 3509 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3510 PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3511 3512 /* call Matlab function in ctx with arguments x and y */ 3513 3514 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 3515 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 3516 ierr = PetscMemcpy(&lxdot,&xdot,sizeof(x));CHKERRQ(ierr); 3517 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 3518 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 3519 prhs[0] = mxCreateDoubleScalar((double)ls); 3520 prhs[1] = mxCreateDoubleScalar((double)time); 3521 prhs[2] = mxCreateDoubleScalar((double)lx); 3522 prhs[3] = mxCreateDoubleScalar((double)lxdot); 3523 prhs[4] = mxCreateDoubleScalar((double)shift); 3524 prhs[5] = mxCreateDoubleScalar((double)lA); 3525 prhs[6] = mxCreateDoubleScalar((double)lB); 3526 prhs[7] = mxCreateString(sctx->funcname); 3527 prhs[8] = sctx->ctx; 3528 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");CHKERRQ(ierr); 3529 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 3530 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 3531 mxDestroyArray(prhs[0]); 3532 mxDestroyArray(prhs[1]); 3533 mxDestroyArray(prhs[2]); 3534 mxDestroyArray(prhs[3]); 3535 mxDestroyArray(prhs[4]); 3536 mxDestroyArray(prhs[5]); 3537 mxDestroyArray(prhs[6]); 3538 mxDestroyArray(prhs[7]); 3539 mxDestroyArray(plhs[0]); 3540 mxDestroyArray(plhs[1]); 3541 PetscFunctionReturn(0); 3542 } 3543 3544 3545 #undef __FUNCT__ 3546 #define __FUNCT__ "TSSetJacobianMatlab" 3547 /* 3548 TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 3549 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 3550 3551 Logically Collective on TS 3552 3553 Input Parameters: 3554 + ts - the TS context 3555 . A,B - Jacobian matrices 3556 . func - function evaluation routine 3557 - ctx - user context 3558 3559 Calling sequence of func: 3560 $ flag = func (TS ts,PetscReal time,Vec x,Vec xdot,Mat A,Mat B,void *ctx); 3561 3562 3563 Level: developer 3564 3565 .keywords: TS, nonlinear, set, function 3566 3567 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 3568 */ 3569 PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx) 3570 { 3571 PetscErrorCode ierr; 3572 TSMatlabContext *sctx; 3573 3574 PetscFunctionBegin; 3575 /* currently sctx is memory bleed */ 3576 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 3577 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 3578 /* 3579 This should work, but it doesn't 3580 sctx->ctx = ctx; 3581 mexMakeArrayPersistent(sctx->ctx); 3582 */ 3583 sctx->ctx = mxDuplicateArray(ctx); 3584 ierr = TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 3585 PetscFunctionReturn(0); 3586 } 3587 3588 #undef __FUNCT__ 3589 #define __FUNCT__ "TSMonitor_Matlab" 3590 /* 3591 TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab(). 3592 3593 Collective on TS 3594 3595 .seealso: TSSetFunction(), TSGetFunction() 3596 @*/ 3597 PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec x, void *ctx) 3598 { 3599 PetscErrorCode ierr; 3600 TSMatlabContext *sctx = (TSMatlabContext *)ctx; 3601 int nlhs = 1,nrhs = 6; 3602 mxArray *plhs[1],*prhs[6]; 3603 long long int lx = 0,ls = 0; 3604 3605 PetscFunctionBegin; 3606 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3607 PetscValidHeaderSpecific(x,VEC_CLASSID,4); 3608 3609 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 3610 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 3611 prhs[0] = mxCreateDoubleScalar((double)ls); 3612 prhs[1] = mxCreateDoubleScalar((double)it); 3613 prhs[2] = mxCreateDoubleScalar((double)time); 3614 prhs[3] = mxCreateDoubleScalar((double)lx); 3615 prhs[4] = mxCreateString(sctx->funcname); 3616 prhs[5] = sctx->ctx; 3617 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");CHKERRQ(ierr); 3618 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 3619 mxDestroyArray(prhs[0]); 3620 mxDestroyArray(prhs[1]); 3621 mxDestroyArray(prhs[2]); 3622 mxDestroyArray(prhs[3]); 3623 mxDestroyArray(prhs[4]); 3624 mxDestroyArray(plhs[0]); 3625 PetscFunctionReturn(0); 3626 } 3627 3628 3629 #undef __FUNCT__ 3630 #define __FUNCT__ "TSMonitorSetMatlab" 3631 /* 3632 TSMonitorSetMatlab - Sets the monitor function from Matlab 3633 3634 Level: developer 3635 3636 .keywords: TS, nonlinear, set, function 3637 3638 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 3639 */ 3640 PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx) 3641 { 3642 PetscErrorCode ierr; 3643 TSMatlabContext *sctx; 3644 3645 PetscFunctionBegin; 3646 /* currently sctx is memory bleed */ 3647 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 3648 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 3649 /* 3650 This should work, but it doesn't 3651 sctx->ctx = ctx; 3652 mexMakeArrayPersistent(sctx->ctx); 3653 */ 3654 sctx->ctx = mxDuplicateArray(ctx); 3655 ierr = TSMonitorSet(ts,TSMonitor_Matlab,sctx,PETSC_NULL);CHKERRQ(ierr); 3656 PetscFunctionReturn(0); 3657 } 3658 #endif 3659