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