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