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 PetscViewerFormat format; 2331 2332 PetscFunctionBegin; 2333 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2334 if (u) PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2335 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 */ 2336 if (!ts->vec_sol || u == ts->vec_sol) { 2337 Vec y; 2338 ierr = VecDuplicate(u,&y);CHKERRQ(ierr); 2339 ierr = TSSetSolution(ts,y);CHKERRQ(ierr); 2340 ierr = VecDestroy(&y);CHKERRQ(ierr); /* grant ownership */ 2341 } 2342 if (u) { 2343 ierr = VecCopy(u,ts->vec_sol);CHKERRQ(ierr); 2344 } 2345 } else { 2346 if (u) { 2347 ierr = TSSetSolution(ts,u);CHKERRQ(ierr); 2348 } 2349 } 2350 ierr = TSSetUp(ts);CHKERRQ(ierr); 2351 /* reset time step and iteration counters */ 2352 ts->steps = 0; 2353 ts->ksp_its = 0; 2354 ts->snes_its = 0; 2355 ts->num_snes_failures = 0; 2356 ts->reject = 0; 2357 ts->reason = TS_CONVERGED_ITERATING; 2358 2359 if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */ 2360 ierr = (*ts->ops->solve)(ts);CHKERRQ(ierr); 2361 ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr); 2362 ts->solvetime = ts->ptime; 2363 } else { 2364 /* steps the requested number of timesteps. */ 2365 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 2366 if (ts->steps >= ts->max_steps) 2367 ts->reason = TS_CONVERGED_ITS; 2368 else if (ts->ptime >= ts->max_time) 2369 ts->reason = TS_CONVERGED_TIME; 2370 while (!ts->reason) { 2371 ierr = TSStep(ts);CHKERRQ(ierr); 2372 ierr = TSPostStep(ts);CHKERRQ(ierr); 2373 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 2374 } 2375 if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) { 2376 ierr = TSInterpolate(ts,ts->max_time,u);CHKERRQ(ierr); 2377 ts->solvetime = ts->max_time; 2378 } else { 2379 if (u) {ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr);} 2380 ts->solvetime = ts->ptime; 2381 } 2382 } 2383 ierr = TSMonitor(ts,-1,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 2384 ierr = PetscOptionsGetViewer(((PetscObject)ts)->comm,((PetscObject)ts)->prefix,"-ts_view",&viewer,&format,&flg);CHKERRQ(ierr); 2385 if (flg && !PetscPreLoadingOn) { 2386 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 2387 ierr = TSView(ts,viewer);CHKERRQ(ierr); 2388 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 2389 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2390 } 2391 PetscFunctionReturn(0); 2392 } 2393 2394 #undef __FUNCT__ 2395 #define __FUNCT__ "TSMonitor" 2396 /*@ 2397 TSMonitor - Runs all user-provided monitor routines set using TSMonitorSet() 2398 2399 Collective on TS 2400 2401 Input Parameters: 2402 + ts - time stepping context obtained from TSCreate() 2403 . step - step number that has just completed 2404 . ptime - model time of the state 2405 - u - state at the current model time 2406 2407 Notes: 2408 TSMonitor() is typically used within the time stepping implementations. 2409 Users might call this function when using the TSStep() interface instead of TSSolve(). 2410 2411 Level: advanced 2412 2413 .keywords: TS, timestep 2414 @*/ 2415 PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u) 2416 { 2417 PetscErrorCode ierr; 2418 PetscInt i,n = ts->numbermonitors; 2419 2420 PetscFunctionBegin; 2421 for (i=0; i<n; i++) { 2422 ierr = (*ts->monitor[i])(ts,step,ptime,u,ts->monitorcontext[i]);CHKERRQ(ierr); 2423 } 2424 PetscFunctionReturn(0); 2425 } 2426 2427 /* ------------------------------------------------------------------------*/ 2428 struct _n_TSMonitorLGCtx { 2429 PetscDrawLG lg; 2430 PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */ 2431 PetscInt ksp_its,snes_its; 2432 }; 2433 2434 2435 #undef __FUNCT__ 2436 #define __FUNCT__ "TSMonitorLGCtxCreate" 2437 /*@C 2438 TSMonitorLGCtxCreate - Creates a line graph context for use with 2439 TS to monitor the solution process graphically in various ways 2440 2441 Collective on TS 2442 2443 Input Parameters: 2444 + host - the X display to open, or null for the local machine 2445 . label - the title to put in the title bar 2446 . x, y - the screen coordinates of the upper left coordinate of the window 2447 . m, n - the screen width and height in pixels 2448 - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time 2449 2450 Output Parameter: 2451 . ctx - the context 2452 2453 Options Database Key: 2454 + -ts_monitor_lg_timestep - automatically sets line graph monitor 2455 . -ts_monitor_lg_solution - 2456 . -ts_monitor_lg_error - 2457 . -ts_monitor_lg_ksp_iterations - 2458 . -ts_monitor_lg_snes_iterations - 2459 - -lg_indicate_data_points <true,false> - indicate the data points (at each time step) on the plot; default is true 2460 2461 Notes: 2462 Use TSMonitorLGCtxDestroy() to destroy. 2463 2464 Level: intermediate 2465 2466 .keywords: TS, monitor, line graph, residual, seealso 2467 2468 .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError() 2469 2470 @*/ 2471 PetscErrorCode TSMonitorLGCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorLGCtx *ctx) 2472 { 2473 PetscDraw win; 2474 PetscErrorCode ierr; 2475 PetscBool flg = PETSC_TRUE; 2476 2477 PetscFunctionBegin; 2478 ierr = PetscNew(struct _n_TSMonitorLGCtx,ctx);CHKERRQ(ierr); 2479 ierr = PetscDrawCreate(comm,host,label,x,y,m,n,&win);CHKERRQ(ierr); 2480 ierr = PetscDrawSetFromOptions(win);CHKERRQ(ierr); 2481 ierr = PetscDrawLGCreate(win,1,&(*ctx)->lg);CHKERRQ(ierr); 2482 ierr = PetscOptionsGetBool(PETSC_NULL,"-lg_indicate_data_points",&flg,PETSC_NULL);CHKERRQ(ierr); 2483 if (flg) { 2484 ierr = PetscDrawLGIndicateDataPoints((*ctx)->lg);CHKERRQ(ierr); 2485 } 2486 ierr = PetscLogObjectParent((*ctx)->lg,win);CHKERRQ(ierr); 2487 (*ctx)->howoften = howoften; 2488 PetscFunctionReturn(0); 2489 } 2490 2491 #undef __FUNCT__ 2492 #define __FUNCT__ "TSMonitorLGTimeStep" 2493 PetscErrorCode TSMonitorLGTimeStep(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 2494 { 2495 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 2496 PetscReal x = ptime,y; 2497 PetscErrorCode ierr; 2498 2499 PetscFunctionBegin; 2500 if (!n) { 2501 PetscDrawAxis axis; 2502 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 2503 ierr = PetscDrawAxisSetLabels(axis,"Timestep as function of time","Time","Time step");CHKERRQ(ierr); 2504 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 2505 } 2506 ierr = TSGetTimeStep(ts,&y);CHKERRQ(ierr); 2507 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 2508 if (((ctx->howoften > 0) && (!(n % ctx->howoften))) || ((ctx->howoften == -1) && (n == -1))){ 2509 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 2510 } 2511 PetscFunctionReturn(0); 2512 } 2513 2514 #undef __FUNCT__ 2515 #define __FUNCT__ "TSMonitorLGCtxDestroy" 2516 /*@C 2517 TSMonitorLGCtxDestroy - Destroys a line graph context that was created 2518 with TSMonitorLGCtxCreate(). 2519 2520 Collective on TSMonitorLGCtx 2521 2522 Input Parameter: 2523 . ctx - the monitor context 2524 2525 Level: intermediate 2526 2527 .keywords: TS, monitor, line graph, destroy 2528 2529 .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep(); 2530 @*/ 2531 PetscErrorCode TSMonitorLGCtxDestroy(TSMonitorLGCtx *ctx) 2532 { 2533 PetscDraw draw; 2534 PetscErrorCode ierr; 2535 2536 PetscFunctionBegin; 2537 ierr = PetscDrawLGGetDraw((*ctx)->lg,&draw);CHKERRQ(ierr); 2538 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 2539 ierr = PetscDrawLGDestroy(&(*ctx)->lg);CHKERRQ(ierr); 2540 ierr = PetscFree(*ctx);CHKERRQ(ierr); 2541 PetscFunctionReturn(0); 2542 } 2543 2544 #undef __FUNCT__ 2545 #define __FUNCT__ "TSGetTime" 2546 /*@ 2547 TSGetTime - Gets the time of the most recently completed step. 2548 2549 Not Collective 2550 2551 Input Parameter: 2552 . ts - the TS context obtained from TSCreate() 2553 2554 Output Parameter: 2555 . t - the current time 2556 2557 Level: beginner 2558 2559 Note: 2560 When called during time step evaluation (e.g. during residual evaluation or via hooks set using TSSetPreStep(), 2561 TSSetPreStage(), or TSSetPostStep()), the time is the time at the start of the step being evaluated. 2562 2563 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 2564 2565 .keywords: TS, get, time 2566 @*/ 2567 PetscErrorCode TSGetTime(TS ts,PetscReal* t) 2568 { 2569 PetscFunctionBegin; 2570 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2571 PetscValidRealPointer(t,2); 2572 *t = ts->ptime; 2573 PetscFunctionReturn(0); 2574 } 2575 2576 #undef __FUNCT__ 2577 #define __FUNCT__ "TSSetTime" 2578 /*@ 2579 TSSetTime - Allows one to reset the time. 2580 2581 Logically Collective on TS 2582 2583 Input Parameters: 2584 + ts - the TS context obtained from TSCreate() 2585 - time - the time 2586 2587 Level: intermediate 2588 2589 .seealso: TSGetTime(), TSSetDuration() 2590 2591 .keywords: TS, set, time 2592 @*/ 2593 PetscErrorCode TSSetTime(TS ts, PetscReal t) 2594 { 2595 PetscFunctionBegin; 2596 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2597 PetscValidLogicalCollectiveReal(ts,t,2); 2598 ts->ptime = t; 2599 PetscFunctionReturn(0); 2600 } 2601 2602 #undef __FUNCT__ 2603 #define __FUNCT__ "TSSetOptionsPrefix" 2604 /*@C 2605 TSSetOptionsPrefix - Sets the prefix used for searching for all 2606 TS options in the database. 2607 2608 Logically Collective on TS 2609 2610 Input Parameter: 2611 + ts - The TS context 2612 - prefix - The prefix to prepend to all option names 2613 2614 Notes: 2615 A hyphen (-) must NOT be given at the beginning of the prefix name. 2616 The first character of all runtime options is AUTOMATICALLY the 2617 hyphen. 2618 2619 Level: advanced 2620 2621 .keywords: TS, set, options, prefix, database 2622 2623 .seealso: TSSetFromOptions() 2624 2625 @*/ 2626 PetscErrorCode TSSetOptionsPrefix(TS ts,const char prefix[]) 2627 { 2628 PetscErrorCode ierr; 2629 SNES snes; 2630 2631 PetscFunctionBegin; 2632 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2633 ierr = PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2634 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2635 ierr = SNESSetOptionsPrefix(snes,prefix);CHKERRQ(ierr); 2636 PetscFunctionReturn(0); 2637 } 2638 2639 2640 #undef __FUNCT__ 2641 #define __FUNCT__ "TSAppendOptionsPrefix" 2642 /*@C 2643 TSAppendOptionsPrefix - Appends to the prefix used for searching for all 2644 TS options in the database. 2645 2646 Logically Collective on TS 2647 2648 Input Parameter: 2649 + ts - The TS context 2650 - prefix - The prefix to prepend to all option names 2651 2652 Notes: 2653 A hyphen (-) must NOT be given at the beginning of the prefix name. 2654 The first character of all runtime options is AUTOMATICALLY the 2655 hyphen. 2656 2657 Level: advanced 2658 2659 .keywords: TS, append, options, prefix, database 2660 2661 .seealso: TSGetOptionsPrefix() 2662 2663 @*/ 2664 PetscErrorCode TSAppendOptionsPrefix(TS ts,const char prefix[]) 2665 { 2666 PetscErrorCode ierr; 2667 SNES snes; 2668 2669 PetscFunctionBegin; 2670 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2671 ierr = PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2672 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2673 ierr = SNESAppendOptionsPrefix(snes,prefix);CHKERRQ(ierr); 2674 PetscFunctionReturn(0); 2675 } 2676 2677 #undef __FUNCT__ 2678 #define __FUNCT__ "TSGetOptionsPrefix" 2679 /*@C 2680 TSGetOptionsPrefix - Sets the prefix used for searching for all 2681 TS options in the database. 2682 2683 Not Collective 2684 2685 Input Parameter: 2686 . ts - The TS context 2687 2688 Output Parameter: 2689 . prefix - A pointer to the prefix string used 2690 2691 Notes: On the fortran side, the user should pass in a string 'prifix' of 2692 sufficient length to hold the prefix. 2693 2694 Level: intermediate 2695 2696 .keywords: TS, get, options, prefix, database 2697 2698 .seealso: TSAppendOptionsPrefix() 2699 @*/ 2700 PetscErrorCode TSGetOptionsPrefix(TS ts,const char *prefix[]) 2701 { 2702 PetscErrorCode ierr; 2703 2704 PetscFunctionBegin; 2705 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2706 PetscValidPointer(prefix,2); 2707 ierr = PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2708 PetscFunctionReturn(0); 2709 } 2710 2711 #undef __FUNCT__ 2712 #define __FUNCT__ "TSGetRHSJacobian" 2713 /*@C 2714 TSGetRHSJacobian - Returns the Jacobian J at the present timestep. 2715 2716 Not Collective, but parallel objects are returned if TS is parallel 2717 2718 Input Parameter: 2719 . ts - The TS context obtained from TSCreate() 2720 2721 Output Parameters: 2722 + J - The Jacobian J of F, where U_t = G(U,t) 2723 . M - The preconditioner matrix, usually the same as J 2724 . func - Function to compute the Jacobian of the RHS 2725 - ctx - User-defined context for Jacobian evaluation routine 2726 2727 Notes: You can pass in PETSC_NULL for any return argument you do not need. 2728 2729 Level: intermediate 2730 2731 .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 2732 2733 .keywords: TS, timestep, get, matrix, Jacobian 2734 @*/ 2735 PetscErrorCode TSGetRHSJacobian(TS ts,Mat *J,Mat *M,TSRHSJacobian *func,void **ctx) 2736 { 2737 PetscErrorCode ierr; 2738 SNES snes; 2739 DM dm; 2740 2741 PetscFunctionBegin; 2742 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2743 ierr = SNESGetJacobian(snes,J,M,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2744 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 2745 ierr = DMTSGetRHSJacobian(dm,func,ctx);CHKERRQ(ierr); 2746 PetscFunctionReturn(0); 2747 } 2748 2749 #undef __FUNCT__ 2750 #define __FUNCT__ "TSGetIJacobian" 2751 /*@C 2752 TSGetIJacobian - Returns the implicit Jacobian at the present timestep. 2753 2754 Not Collective, but parallel objects are returned if TS is parallel 2755 2756 Input Parameter: 2757 . ts - The TS context obtained from TSCreate() 2758 2759 Output Parameters: 2760 + A - The Jacobian of F(t,U,U_t) 2761 . B - The preconditioner matrix, often the same as A 2762 . f - The function to compute the matrices 2763 - ctx - User-defined context for Jacobian evaluation routine 2764 2765 Notes: You can pass in PETSC_NULL for any return argument you do not need. 2766 2767 Level: advanced 2768 2769 .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 2770 2771 .keywords: TS, timestep, get, matrix, Jacobian 2772 @*/ 2773 PetscErrorCode TSGetIJacobian(TS ts,Mat *A,Mat *B,TSIJacobian *f,void **ctx) 2774 { 2775 PetscErrorCode ierr; 2776 SNES snes; 2777 DM dm; 2778 2779 PetscFunctionBegin; 2780 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2781 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 2782 ierr = SNESGetJacobian(snes,A,B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2783 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 2784 ierr = DMTSGetIJacobian(dm,f,ctx);CHKERRQ(ierr); 2785 PetscFunctionReturn(0); 2786 } 2787 2788 struct _n_TSMonitorDrawCtx { 2789 PetscViewer viewer; 2790 Vec initialsolution; 2791 PetscBool showinitial; 2792 PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */ 2793 }; 2794 2795 #undef __FUNCT__ 2796 #define __FUNCT__ "TSMonitorDrawSolution" 2797 /*@C 2798 TSMonitorDrawSolution - Monitors progress of the TS solvers by calling 2799 VecView() for the solution at each timestep 2800 2801 Collective on TS 2802 2803 Input Parameters: 2804 + ts - the TS context 2805 . step - current time-step 2806 . ptime - current time 2807 - dummy - either a viewer or PETSC_NULL 2808 2809 Options Database: 2810 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 2811 2812 Notes: the initial solution and current solution are not displayed with a common axis scaling so generally the option -ts_monitor_draw_solution_initial 2813 will look bad 2814 2815 Level: intermediate 2816 2817 .keywords: TS, vector, monitor, view 2818 2819 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 2820 @*/ 2821 PetscErrorCode TSMonitorDrawSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 2822 { 2823 PetscErrorCode ierr; 2824 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy; 2825 2826 PetscFunctionBegin; 2827 if (!step && ictx->showinitial) { 2828 if (!ictx->initialsolution) { 2829 ierr = VecDuplicate(u,&ictx->initialsolution);CHKERRQ(ierr); 2830 } 2831 ierr = VecCopy(u,ictx->initialsolution);CHKERRQ(ierr); 2832 } 2833 if (!(((ictx->howoften > 0) && (!(step % ictx->howoften)) && (step > -1)) || ((ictx->howoften == -1) && (step == -1)))) PetscFunctionReturn(0); 2834 2835 if (ictx->showinitial) { 2836 PetscReal pause; 2837 ierr = PetscViewerDrawGetPause(ictx->viewer,&pause);CHKERRQ(ierr); 2838 ierr = PetscViewerDrawSetPause(ictx->viewer,0.0);CHKERRQ(ierr); 2839 ierr = VecView(ictx->initialsolution,ictx->viewer);CHKERRQ(ierr); 2840 ierr = PetscViewerDrawSetPause(ictx->viewer,pause);CHKERRQ(ierr); 2841 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);CHKERRQ(ierr); 2842 } 2843 ierr = VecView(u,ictx->viewer);CHKERRQ(ierr); 2844 if (ictx->showinitial) { 2845 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);CHKERRQ(ierr); 2846 } 2847 PetscFunctionReturn(0); 2848 } 2849 2850 2851 #undef __FUNCT__ 2852 #define __FUNCT__ "TSMonitorDrawCtxDestroy" 2853 /*@C 2854 TSMonitorDrawCtxDestroy - Destroys the monitor context for TSMonitorDrawSolution() 2855 2856 Collective on TS 2857 2858 Input Parameters: 2859 . ctx - the monitor context 2860 2861 Level: intermediate 2862 2863 .keywords: TS, vector, monitor, view 2864 2865 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawSolution(), TSMonitorDrawError() 2866 @*/ 2867 PetscErrorCode TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx) 2868 { 2869 PetscErrorCode ierr; 2870 2871 PetscFunctionBegin; 2872 ierr = PetscViewerDestroy(&(*ictx)->viewer);CHKERRQ(ierr); 2873 ierr = VecDestroy(&(*ictx)->initialsolution);CHKERRQ(ierr); 2874 ierr = PetscFree(*ictx);CHKERRQ(ierr); 2875 PetscFunctionReturn(0); 2876 } 2877 2878 #undef __FUNCT__ 2879 #define __FUNCT__ "TSMonitorDrawCtxCreate" 2880 /*@C 2881 TSMonitorDrawCtxCreate - Creates the monitor context for TSMonitorDrawCtx 2882 2883 Collective on TS 2884 2885 Input Parameter: 2886 . ts - time-step context 2887 2888 Output Patameter: 2889 . ctx - the monitor context 2890 2891 Options Database: 2892 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 2893 2894 Level: intermediate 2895 2896 .keywords: TS, vector, monitor, view 2897 2898 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawCtx() 2899 @*/ 2900 PetscErrorCode TSMonitorDrawCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorDrawCtx *ctx) 2901 { 2902 PetscErrorCode ierr; 2903 2904 PetscFunctionBegin; 2905 ierr = PetscNew(struct _n_TSMonitorDrawCtx,ctx);CHKERRQ(ierr); 2906 ierr = PetscViewerDrawOpen(comm,host,label,x,y,m,n,&(*ctx)->viewer);CHKERRQ(ierr); 2907 (*ctx)->showinitial = PETSC_FALSE; 2908 (*ctx)->howoften = howoften; 2909 ierr = PetscOptionsGetBool(PETSC_NULL,"-ts_monitor_draw_solution_initial",&(*ctx)->showinitial,PETSC_NULL);CHKERRQ(ierr); 2910 PetscFunctionReturn(0); 2911 } 2912 2913 #undef __FUNCT__ 2914 #define __FUNCT__ "TSMonitorDrawError" 2915 /*@C 2916 TSMonitorDrawError - Monitors progress of the TS solvers by calling 2917 VecView() for the error at each timestep 2918 2919 Collective on TS 2920 2921 Input Parameters: 2922 + ts - the TS context 2923 . step - current time-step 2924 . ptime - current time 2925 - dummy - either a viewer or PETSC_NULL 2926 2927 Level: intermediate 2928 2929 .keywords: TS, vector, monitor, view 2930 2931 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 2932 @*/ 2933 PetscErrorCode TSMonitorDrawError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 2934 { 2935 PetscErrorCode ierr; 2936 TSMonitorDrawCtx ctx = (TSMonitorDrawCtx)dummy; 2937 PetscViewer viewer = ctx->viewer; 2938 Vec work; 2939 2940 PetscFunctionBegin; 2941 if (!(((ctx->howoften > 0) && (!(step % ctx->howoften)) && (step > -1)) || ((ctx->howoften == -1) && (step == -1)))) PetscFunctionReturn(0); 2942 ierr = VecDuplicate(u,&work);CHKERRQ(ierr); 2943 ierr = TSComputeSolutionFunction(ts,ptime,work);CHKERRQ(ierr); 2944 ierr = VecAXPY(work,-1.0,u);CHKERRQ(ierr); 2945 ierr = VecView(work,viewer);CHKERRQ(ierr); 2946 ierr = VecDestroy(&work);CHKERRQ(ierr); 2947 PetscFunctionReturn(0); 2948 } 2949 2950 #include <petsc-private/dmimpl.h> 2951 #undef __FUNCT__ 2952 #define __FUNCT__ "TSSetDM" 2953 /*@ 2954 TSSetDM - Sets the DM that may be used by some preconditioners 2955 2956 Logically Collective on TS and DM 2957 2958 Input Parameters: 2959 + ts - the preconditioner context 2960 - dm - the dm 2961 2962 Level: intermediate 2963 2964 2965 .seealso: TSGetDM(), SNESSetDM(), SNESGetDM() 2966 @*/ 2967 PetscErrorCode TSSetDM(TS ts,DM dm) 2968 { 2969 PetscErrorCode ierr; 2970 SNES snes; 2971 DMTS tsdm; 2972 2973 PetscFunctionBegin; 2974 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2975 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 2976 if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */ 2977 if (ts->dm->dmts && !dm->dmts) { 2978 ierr = DMCopyDMTS(ts->dm,dm);CHKERRQ(ierr); 2979 ierr = DMGetDMTS(ts->dm,&tsdm);CHKERRQ(ierr); 2980 if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */ 2981 tsdm->originaldm = dm; 2982 } 2983 } 2984 ierr = DMDestroy(&ts->dm);CHKERRQ(ierr); 2985 } 2986 ts->dm = dm; 2987 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2988 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 2989 PetscFunctionReturn(0); 2990 } 2991 2992 #undef __FUNCT__ 2993 #define __FUNCT__ "TSGetDM" 2994 /*@ 2995 TSGetDM - Gets the DM that may be used by some preconditioners 2996 2997 Not Collective 2998 2999 Input Parameter: 3000 . ts - the preconditioner context 3001 3002 Output Parameter: 3003 . dm - the dm 3004 3005 Level: intermediate 3006 3007 3008 .seealso: TSSetDM(), SNESSetDM(), SNESGetDM() 3009 @*/ 3010 PetscErrorCode TSGetDM(TS ts,DM *dm) 3011 { 3012 PetscErrorCode ierr; 3013 3014 PetscFunctionBegin; 3015 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3016 if (!ts->dm) { 3017 ierr = DMShellCreate(((PetscObject)ts)->comm,&ts->dm);CHKERRQ(ierr); 3018 if (ts->snes) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 3019 } 3020 *dm = ts->dm; 3021 PetscFunctionReturn(0); 3022 } 3023 3024 #undef __FUNCT__ 3025 #define __FUNCT__ "SNESTSFormFunction" 3026 /*@ 3027 SNESTSFormFunction - Function to evaluate nonlinear residual 3028 3029 Logically Collective on SNES 3030 3031 Input Parameter: 3032 + snes - nonlinear solver 3033 . U - the current state at which to evaluate the residual 3034 - ctx - user context, must be a TS 3035 3036 Output Parameter: 3037 . F - the nonlinear residual 3038 3039 Notes: 3040 This function is not normally called by users and is automatically registered with the SNES used by TS. 3041 It is most frequently passed to MatFDColoringSetFunction(). 3042 3043 Level: advanced 3044 3045 .seealso: SNESSetFunction(), MatFDColoringSetFunction() 3046 @*/ 3047 PetscErrorCode SNESTSFormFunction(SNES snes,Vec U,Vec F,void *ctx) 3048 { 3049 TS ts = (TS)ctx; 3050 PetscErrorCode ierr; 3051 3052 PetscFunctionBegin; 3053 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3054 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 3055 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 3056 PetscValidHeaderSpecific(ts,TS_CLASSID,4); 3057 ierr = (ts->ops->snesfunction)(snes,U,F,ts);CHKERRQ(ierr); 3058 PetscFunctionReturn(0); 3059 } 3060 3061 #undef __FUNCT__ 3062 #define __FUNCT__ "SNESTSFormJacobian" 3063 /*@ 3064 SNESTSFormJacobian - Function to evaluate the Jacobian 3065 3066 Collective on SNES 3067 3068 Input Parameter: 3069 + snes - nonlinear solver 3070 . U - the current state at which to evaluate the residual 3071 - ctx - user context, must be a TS 3072 3073 Output Parameter: 3074 + A - the Jacobian 3075 . B - the preconditioning matrix (may be the same as A) 3076 - flag - indicates any structure change in the matrix 3077 3078 Notes: 3079 This function is not normally called by users and is automatically registered with the SNES used by TS. 3080 3081 Level: developer 3082 3083 .seealso: SNESSetJacobian() 3084 @*/ 3085 PetscErrorCode SNESTSFormJacobian(SNES snes,Vec U,Mat *A,Mat *B,MatStructure *flag,void *ctx) 3086 { 3087 TS ts = (TS)ctx; 3088 PetscErrorCode ierr; 3089 3090 PetscFunctionBegin; 3091 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3092 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 3093 PetscValidPointer(A,3); 3094 PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 3095 PetscValidPointer(B,4); 3096 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); 3097 PetscValidPointer(flag,5); 3098 PetscValidHeaderSpecific(ts,TS_CLASSID,6); 3099 ierr = (ts->ops->snesjacobian)(snes,U,A,B,flag,ts);CHKERRQ(ierr); 3100 PetscFunctionReturn(0); 3101 } 3102 3103 #undef __FUNCT__ 3104 #define __FUNCT__ "TSComputeRHSFunctionLinear" 3105 /*@C 3106 TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems only 3107 3108 Collective on TS 3109 3110 Input Arguments: 3111 + ts - time stepping context 3112 . t - time at which to evaluate 3113 . U - state at which to evaluate 3114 - ctx - context 3115 3116 Output Arguments: 3117 . F - right hand side 3118 3119 Level: intermediate 3120 3121 Notes: 3122 This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems. 3123 The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian(). 3124 3125 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant() 3126 @*/ 3127 PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec U,Vec F,void *ctx) 3128 { 3129 PetscErrorCode ierr; 3130 Mat Arhs,Brhs; 3131 MatStructure flg2; 3132 3133 PetscFunctionBegin; 3134 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 3135 ierr = TSComputeRHSJacobian(ts,t,U,&Arhs,&Brhs,&flg2);CHKERRQ(ierr); 3136 ierr = MatMult(Arhs,U,F);CHKERRQ(ierr); 3137 PetscFunctionReturn(0); 3138 } 3139 3140 #undef __FUNCT__ 3141 #define __FUNCT__ "TSComputeRHSJacobianConstant" 3142 /*@C 3143 TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent. 3144 3145 Collective on TS 3146 3147 Input Arguments: 3148 + ts - time stepping context 3149 . t - time at which to evaluate 3150 . U - state at which to evaluate 3151 - ctx - context 3152 3153 Output Arguments: 3154 + A - pointer to operator 3155 . B - pointer to preconditioning matrix 3156 - flg - matrix structure flag 3157 3158 Level: intermediate 3159 3160 Notes: 3161 This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems. 3162 3163 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear() 3164 @*/ 3165 PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec U,Mat *A,Mat *B,MatStructure *flg,void *ctx) 3166 { 3167 PetscFunctionBegin; 3168 *flg = SAME_PRECONDITIONER; 3169 PetscFunctionReturn(0); 3170 } 3171 3172 #undef __FUNCT__ 3173 #define __FUNCT__ "TSComputeIFunctionLinear" 3174 /*@C 3175 TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only 3176 3177 Collective on TS 3178 3179 Input Arguments: 3180 + ts - time stepping context 3181 . t - time at which to evaluate 3182 . U - state at which to evaluate 3183 . Udot - time derivative of state vector 3184 - ctx - context 3185 3186 Output Arguments: 3187 . F - left hand side 3188 3189 Level: intermediate 3190 3191 Notes: 3192 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 3193 user is required to write their own TSComputeIFunction. 3194 This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems. 3195 The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian(). 3196 3197 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant() 3198 @*/ 3199 PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx) 3200 { 3201 PetscErrorCode ierr; 3202 Mat A,B; 3203 MatStructure flg2; 3204 3205 PetscFunctionBegin; 3206 ierr = TSGetIJacobian(ts,&A,&B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 3207 ierr = TSComputeIJacobian(ts,t,U,Udot,1.0,&A,&B,&flg2,PETSC_TRUE);CHKERRQ(ierr); 3208 ierr = MatMult(A,Udot,F);CHKERRQ(ierr); 3209 PetscFunctionReturn(0); 3210 } 3211 3212 #undef __FUNCT__ 3213 #define __FUNCT__ "TSComputeIJacobianConstant" 3214 /*@C 3215 TSComputeIJacobianConstant - Reuses a Jacobian that is time-independent. 3216 3217 Collective on TS 3218 3219 Input Arguments: 3220 + ts - time stepping context 3221 . t - time at which to evaluate 3222 . U - state at which to evaluate 3223 . Udot - time derivative of state vector 3224 . shift - shift to apply 3225 - ctx - context 3226 3227 Output Arguments: 3228 + A - pointer to operator 3229 . B - pointer to preconditioning matrix 3230 - flg - matrix structure flag 3231 3232 Level: intermediate 3233 3234 Notes: 3235 This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems. 3236 3237 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear() 3238 @*/ 3239 PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,void *ctx) 3240 { 3241 PetscFunctionBegin; 3242 *flg = SAME_PRECONDITIONER; 3243 PetscFunctionReturn(0); 3244 } 3245 3246 #undef __FUNCT__ 3247 #define __FUNCT__ "TSGetConvergedReason" 3248 /*@ 3249 TSGetConvergedReason - Gets the reason the TS iteration was stopped. 3250 3251 Not Collective 3252 3253 Input Parameter: 3254 . ts - the TS context 3255 3256 Output Parameter: 3257 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 3258 manual pages for the individual convergence tests for complete lists 3259 3260 Level: beginner 3261 3262 Notes: 3263 Can only be called after the call to TSSolve() is complete. 3264 3265 .keywords: TS, nonlinear, set, convergence, test 3266 3267 .seealso: TSSetConvergenceTest(), TSConvergedReason 3268 @*/ 3269 PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason) 3270 { 3271 PetscFunctionBegin; 3272 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3273 PetscValidPointer(reason,2); 3274 *reason = ts->reason; 3275 PetscFunctionReturn(0); 3276 } 3277 3278 #undef __FUNCT__ 3279 #define __FUNCT__ "TSSetConvergedReason" 3280 /*@ 3281 TSSetConvergedReason - Sets the reason for handling the convergence of TSSolve. 3282 3283 Not Collective 3284 3285 Input Parameter: 3286 + ts - the TS context 3287 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 3288 manual pages for the individual convergence tests for complete lists 3289 3290 Level: advanced 3291 3292 Notes: 3293 Can only be called during TSSolve() is active. 3294 3295 .keywords: TS, nonlinear, set, convergence, test 3296 3297 .seealso: TSConvergedReason 3298 @*/ 3299 PetscErrorCode TSSetConvergedReason(TS ts,TSConvergedReason reason) 3300 { 3301 PetscFunctionBegin; 3302 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3303 ts->reason = reason; 3304 PetscFunctionReturn(0); 3305 } 3306 3307 #undef __FUNCT__ 3308 #define __FUNCT__ "TSGetSolveTime" 3309 /*@ 3310 TSGetSolveTime - Gets the time after a call to TSSolve() 3311 3312 Not Collective 3313 3314 Input Parameter: 3315 . ts - the TS context 3316 3317 Output Parameter: 3318 . ftime - the final time. This time should correspond to the final time set with TSSetDuration() 3319 3320 Level: beginner 3321 3322 Notes: 3323 Can only be called after the call to TSSolve() is complete. 3324 3325 .keywords: TS, nonlinear, set, convergence, test 3326 3327 .seealso: TSSetConvergenceTest(), TSConvergedReason 3328 @*/ 3329 PetscErrorCode TSGetSolveTime(TS ts,PetscReal *ftime) 3330 { 3331 PetscFunctionBegin; 3332 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3333 PetscValidPointer(ftime,2); 3334 *ftime = ts->solvetime; 3335 PetscFunctionReturn(0); 3336 } 3337 3338 #undef __FUNCT__ 3339 #define __FUNCT__ "TSGetSNESIterations" 3340 /*@ 3341 TSGetSNESIterations - Gets the total number of nonlinear iterations 3342 used by the time integrator. 3343 3344 Not Collective 3345 3346 Input Parameter: 3347 . ts - TS context 3348 3349 Output Parameter: 3350 . nits - number of nonlinear iterations 3351 3352 Notes: 3353 This counter is reset to zero for each successive call to TSSolve(). 3354 3355 Level: intermediate 3356 3357 .keywords: TS, get, number, nonlinear, iterations 3358 3359 .seealso: TSGetKSPIterations() 3360 @*/ 3361 PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits) 3362 { 3363 PetscFunctionBegin; 3364 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3365 PetscValidIntPointer(nits,2); 3366 *nits = ts->snes_its; 3367 PetscFunctionReturn(0); 3368 } 3369 3370 #undef __FUNCT__ 3371 #define __FUNCT__ "TSGetKSPIterations" 3372 /*@ 3373 TSGetKSPIterations - Gets the total number of linear iterations 3374 used by the time integrator. 3375 3376 Not Collective 3377 3378 Input Parameter: 3379 . ts - TS context 3380 3381 Output Parameter: 3382 . lits - number of linear iterations 3383 3384 Notes: 3385 This counter is reset to zero for each successive call to TSSolve(). 3386 3387 Level: intermediate 3388 3389 .keywords: TS, get, number, linear, iterations 3390 3391 .seealso: TSGetSNESIterations(), SNESGetKSPIterations() 3392 @*/ 3393 PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits) 3394 { 3395 PetscFunctionBegin; 3396 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3397 PetscValidIntPointer(lits,2); 3398 *lits = ts->ksp_its; 3399 PetscFunctionReturn(0); 3400 } 3401 3402 #undef __FUNCT__ 3403 #define __FUNCT__ "TSGetStepRejections" 3404 /*@ 3405 TSGetStepRejections - Gets the total number of rejected steps. 3406 3407 Not Collective 3408 3409 Input Parameter: 3410 . ts - TS context 3411 3412 Output Parameter: 3413 . rejects - number of steps rejected 3414 3415 Notes: 3416 This counter is reset to zero for each successive call to TSSolve(). 3417 3418 Level: intermediate 3419 3420 .keywords: TS, get, number 3421 3422 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails() 3423 @*/ 3424 PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects) 3425 { 3426 PetscFunctionBegin; 3427 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3428 PetscValidIntPointer(rejects,2); 3429 *rejects = ts->reject; 3430 PetscFunctionReturn(0); 3431 } 3432 3433 #undef __FUNCT__ 3434 #define __FUNCT__ "TSGetSNESFailures" 3435 /*@ 3436 TSGetSNESFailures - Gets the total number of failed SNES solves 3437 3438 Not Collective 3439 3440 Input Parameter: 3441 . ts - TS context 3442 3443 Output Parameter: 3444 . fails - number of failed nonlinear solves 3445 3446 Notes: 3447 This counter is reset to zero for each successive call to TSSolve(). 3448 3449 Level: intermediate 3450 3451 .keywords: TS, get, number 3452 3453 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures() 3454 @*/ 3455 PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails) 3456 { 3457 PetscFunctionBegin; 3458 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3459 PetscValidIntPointer(fails,2); 3460 *fails = ts->num_snes_failures; 3461 PetscFunctionReturn(0); 3462 } 3463 3464 #undef __FUNCT__ 3465 #define __FUNCT__ "TSSetMaxStepRejections" 3466 /*@ 3467 TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails 3468 3469 Not Collective 3470 3471 Input Parameter: 3472 + ts - TS context 3473 - rejects - maximum number of rejected steps, pass -1 for unlimited 3474 3475 Notes: 3476 The counter is reset to zero for each step 3477 3478 Options Database Key: 3479 . -ts_max_reject - Maximum number of step rejections before a step fails 3480 3481 Level: intermediate 3482 3483 .keywords: TS, set, maximum, number 3484 3485 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 3486 @*/ 3487 PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects) 3488 { 3489 PetscFunctionBegin; 3490 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3491 ts->max_reject = rejects; 3492 PetscFunctionReturn(0); 3493 } 3494 3495 #undef __FUNCT__ 3496 #define __FUNCT__ "TSSetMaxSNESFailures" 3497 /*@ 3498 TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves 3499 3500 Not Collective 3501 3502 Input Parameter: 3503 + ts - TS context 3504 - fails - maximum number of failed nonlinear solves, pass -1 for unlimited 3505 3506 Notes: 3507 The counter is reset to zero for each successive call to TSSolve(). 3508 3509 Options Database Key: 3510 . -ts_max_snes_failures - Maximum number of nonlinear solve failures 3511 3512 Level: intermediate 3513 3514 .keywords: TS, set, maximum, number 3515 3516 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason() 3517 @*/ 3518 PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails) 3519 { 3520 PetscFunctionBegin; 3521 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3522 ts->max_snes_failures = fails; 3523 PetscFunctionReturn(0); 3524 } 3525 3526 #undef __FUNCT__ 3527 #define __FUNCT__ "TSSetErrorIfStepFails()" 3528 /*@ 3529 TSSetErrorIfStepFails - Error if no step succeeds 3530 3531 Not Collective 3532 3533 Input Parameter: 3534 + ts - TS context 3535 - err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure 3536 3537 Options Database Key: 3538 . -ts_error_if_step_fails - Error if no step succeeds 3539 3540 Level: intermediate 3541 3542 .keywords: TS, set, error 3543 3544 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 3545 @*/ 3546 PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err) 3547 { 3548 PetscFunctionBegin; 3549 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3550 ts->errorifstepfailed = err; 3551 PetscFunctionReturn(0); 3552 } 3553 3554 #undef __FUNCT__ 3555 #define __FUNCT__ "TSMonitorSolutionBinary" 3556 /*@C 3557 TSMonitorSolutionBinary - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file 3558 3559 Collective on TS 3560 3561 Input Parameters: 3562 + ts - the TS context 3563 . step - current time-step 3564 . ptime - current time 3565 . u - current state 3566 - viewer - binary viewer 3567 3568 Level: intermediate 3569 3570 .keywords: TS, vector, monitor, view 3571 3572 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3573 @*/ 3574 PetscErrorCode TSMonitorSolutionBinary(TS ts,PetscInt step,PetscReal ptime,Vec u,void *viewer) 3575 { 3576 PetscErrorCode ierr; 3577 PetscViewer v = (PetscViewer)viewer; 3578 3579 PetscFunctionBegin; 3580 ierr = VecView(u,v);CHKERRQ(ierr); 3581 PetscFunctionReturn(0); 3582 } 3583 3584 #undef __FUNCT__ 3585 #define __FUNCT__ "TSMonitorSolutionVTK" 3586 /*@C 3587 TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep. 3588 3589 Collective on TS 3590 3591 Input Parameters: 3592 + ts - the TS context 3593 . step - current time-step 3594 . ptime - current time 3595 . u - current state 3596 - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 3597 3598 Level: intermediate 3599 3600 Notes: 3601 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. 3602 These are named according to the file name template. 3603 3604 This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy(). 3605 3606 .keywords: TS, vector, monitor, view 3607 3608 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3609 @*/ 3610 PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec u,void *filenametemplate) 3611 { 3612 PetscErrorCode ierr; 3613 char filename[PETSC_MAX_PATH_LEN]; 3614 PetscViewer viewer; 3615 3616 PetscFunctionBegin; 3617 ierr = PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);CHKERRQ(ierr); 3618 ierr = PetscViewerVTKOpen(((PetscObject)ts)->comm,filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 3619 ierr = VecView(u,viewer);CHKERRQ(ierr); 3620 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3621 PetscFunctionReturn(0); 3622 } 3623 3624 #undef __FUNCT__ 3625 #define __FUNCT__ "TSMonitorSolutionVTKDestroy" 3626 /*@C 3627 TSMonitorSolutionVTKDestroy - Destroy context for monitoring 3628 3629 Collective on TS 3630 3631 Input Parameters: 3632 . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 3633 3634 Level: intermediate 3635 3636 Note: 3637 This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK(). 3638 3639 .keywords: TS, vector, monitor, view 3640 3641 .seealso: TSMonitorSet(), TSMonitorSolutionVTK() 3642 @*/ 3643 PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate) 3644 { 3645 PetscErrorCode ierr; 3646 3647 PetscFunctionBegin; 3648 ierr = PetscFree(*(char**)filenametemplate);CHKERRQ(ierr); 3649 PetscFunctionReturn(0); 3650 } 3651 3652 #undef __FUNCT__ 3653 #define __FUNCT__ "TSGetTSAdapt" 3654 /*@ 3655 TSGetTSAdapt - Get the adaptive controller context for the current method 3656 3657 Collective on TS if controller has not been created yet 3658 3659 Input Arguments: 3660 . ts - time stepping context 3661 3662 Output Arguments: 3663 . adapt - adaptive controller 3664 3665 Level: intermediate 3666 3667 .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose() 3668 @*/ 3669 PetscErrorCode TSGetTSAdapt(TS ts,TSAdapt *adapt) 3670 { 3671 PetscErrorCode ierr; 3672 3673 PetscFunctionBegin; 3674 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3675 PetscValidPointer(adapt,2); 3676 if (!ts->adapt) { 3677 ierr = TSAdaptCreate(((PetscObject)ts)->comm,&ts->adapt);CHKERRQ(ierr); 3678 ierr = PetscLogObjectParent(ts,ts->adapt);CHKERRQ(ierr); 3679 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);CHKERRQ(ierr); 3680 } 3681 *adapt = ts->adapt; 3682 PetscFunctionReturn(0); 3683 } 3684 3685 #undef __FUNCT__ 3686 #define __FUNCT__ "TSSetTolerances" 3687 /*@ 3688 TSSetTolerances - Set tolerances for local truncation error when using adaptive controller 3689 3690 Logically Collective 3691 3692 Input Arguments: 3693 + ts - time integration context 3694 . atol - scalar absolute tolerances, PETSC_DECIDE to leave current value 3695 . vatol - vector of absolute tolerances or PETSC_NULL, used in preference to atol if present 3696 . rtol - scalar relative tolerances, PETSC_DECIDE to leave current value 3697 - vrtol - vector of relative tolerances or PETSC_NULL, used in preference to atol if present 3698 3699 Level: beginner 3700 3701 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSGetTolerances() 3702 @*/ 3703 PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol) 3704 { 3705 PetscErrorCode ierr; 3706 3707 PetscFunctionBegin; 3708 if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol; 3709 if (vatol) { 3710 ierr = PetscObjectReference((PetscObject)vatol);CHKERRQ(ierr); 3711 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 3712 ts->vatol = vatol; 3713 } 3714 if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol; 3715 if (vrtol) { 3716 ierr = PetscObjectReference((PetscObject)vrtol);CHKERRQ(ierr); 3717 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 3718 ts->vrtol = vrtol; 3719 } 3720 PetscFunctionReturn(0); 3721 } 3722 3723 #undef __FUNCT__ 3724 #define __FUNCT__ "TSGetTolerances" 3725 /*@ 3726 TSGetTolerances - Get tolerances for local truncation error when using adaptive controller 3727 3728 Logically Collective 3729 3730 Input Arguments: 3731 . ts - time integration context 3732 3733 Output Arguments: 3734 + atol - scalar absolute tolerances, PETSC_NULL to ignore 3735 . vatol - vector of absolute tolerances, PETSC_NULL to ignore 3736 . rtol - scalar relative tolerances, PETSC_NULL to ignore 3737 - vrtol - vector of relative tolerances, PETSC_NULL to ignore 3738 3739 Level: beginner 3740 3741 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSSetTolerances() 3742 @*/ 3743 PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol) 3744 { 3745 PetscFunctionBegin; 3746 if (atol) *atol = ts->atol; 3747 if (vatol) *vatol = ts->vatol; 3748 if (rtol) *rtol = ts->rtol; 3749 if (vrtol) *vrtol = ts->vrtol; 3750 PetscFunctionReturn(0); 3751 } 3752 3753 #undef __FUNCT__ 3754 #define __FUNCT__ "TSErrorNormWRMS" 3755 /*@ 3756 TSErrorNormWRMS - compute a weighted norm of the difference between a vector and the current state 3757 3758 Collective on TS 3759 3760 Input Arguments: 3761 + ts - time stepping context 3762 - Y - state vector to be compared to ts->vec_sol 3763 3764 Output Arguments: 3765 . norm - weighted norm, a value of 1.0 is considered small 3766 3767 Level: developer 3768 3769 .seealso: TSSetTolerances() 3770 @*/ 3771 PetscErrorCode TSErrorNormWRMS(TS ts,Vec Y,PetscReal *norm) 3772 { 3773 PetscErrorCode ierr; 3774 PetscInt i,n,N; 3775 const PetscScalar *u,*y; 3776 Vec U; 3777 PetscReal sum,gsum; 3778 3779 PetscFunctionBegin; 3780 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3781 PetscValidHeaderSpecific(Y,VEC_CLASSID,2); 3782 PetscValidPointer(norm,3); 3783 U = ts->vec_sol; 3784 PetscCheckSameTypeAndComm(U,1,Y,2); 3785 if (U == Y) SETERRQ(((PetscObject)U)->comm,PETSC_ERR_ARG_IDN,"Y cannot be the TS solution vector"); 3786 3787 ierr = VecGetSize(U,&N);CHKERRQ(ierr); 3788 ierr = VecGetLocalSize(U,&n);CHKERRQ(ierr); 3789 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 3790 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 3791 sum = 0.; 3792 if (ts->vatol && ts->vrtol) { 3793 const PetscScalar *atol,*rtol; 3794 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3795 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3796 for (i=0; i<n; i++) { 3797 PetscReal tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 3798 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 3799 } 3800 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3801 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3802 } else if (ts->vatol) { /* vector atol, scalar rtol */ 3803 const PetscScalar *atol; 3804 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3805 for (i=0; i<n; i++) { 3806 PetscReal tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 3807 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 3808 } 3809 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3810 } else if (ts->vrtol) { /* scalar atol, vector rtol */ 3811 const PetscScalar *rtol; 3812 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3813 for (i=0; i<n; i++) { 3814 PetscReal tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 3815 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 3816 } 3817 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3818 } else { /* scalar atol, scalar rtol */ 3819 for (i=0; i<n; i++) { 3820 PetscReal tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 3821 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 3822 } 3823 } 3824 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 3825 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 3826 3827 ierr = MPI_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,((PetscObject)ts)->comm);CHKERRQ(ierr); 3828 *norm = PetscSqrtReal(gsum / N); 3829 if (PetscIsInfOrNanScalar(*norm)) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 3830 PetscFunctionReturn(0); 3831 } 3832 3833 #undef __FUNCT__ 3834 #define __FUNCT__ "TSSetCFLTimeLocal" 3835 /*@ 3836 TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler 3837 3838 Logically Collective on TS 3839 3840 Input Arguments: 3841 + ts - time stepping context 3842 - cfltime - maximum stable time step if using forward Euler (value can be different on each process) 3843 3844 Note: 3845 After calling this function, the global CFL time can be obtained by calling TSGetCFLTime() 3846 3847 Level: intermediate 3848 3849 .seealso: TSGetCFLTime(), TSADAPTCFL 3850 @*/ 3851 PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime) 3852 { 3853 PetscFunctionBegin; 3854 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3855 ts->cfltime_local = cfltime; 3856 ts->cfltime = -1.; 3857 PetscFunctionReturn(0); 3858 } 3859 3860 #undef __FUNCT__ 3861 #define __FUNCT__ "TSGetCFLTime" 3862 /*@ 3863 TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler 3864 3865 Collective on TS 3866 3867 Input Arguments: 3868 . ts - time stepping context 3869 3870 Output Arguments: 3871 . cfltime - maximum stable time step for forward Euler 3872 3873 Level: advanced 3874 3875 .seealso: TSSetCFLTimeLocal() 3876 @*/ 3877 PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime) 3878 { 3879 PetscErrorCode ierr; 3880 3881 PetscFunctionBegin; 3882 if (ts->cfltime < 0) { 3883 ierr = MPI_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,((PetscObject)ts)->comm);CHKERRQ(ierr); 3884 } 3885 *cfltime = ts->cfltime; 3886 PetscFunctionReturn(0); 3887 } 3888 3889 #undef __FUNCT__ 3890 #define __FUNCT__ "TSVISetVariableBounds" 3891 /*@ 3892 TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu 3893 3894 Input Parameters: 3895 . ts - the TS context. 3896 . xl - lower bound. 3897 . xu - upper bound. 3898 3899 Notes: 3900 If this routine is not called then the lower and upper bounds are set to 3901 SNES_VI_NINF and SNES_VI_INF respectively during SNESSetUp(). 3902 3903 Level: advanced 3904 3905 @*/ 3906 PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu) 3907 { 3908 PetscErrorCode ierr; 3909 SNES snes; 3910 3911 PetscFunctionBegin; 3912 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3913 ierr = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(ierr); 3914 PetscFunctionReturn(0); 3915 } 3916 3917 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3918 #include <mex.h> 3919 3920 typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext; 3921 3922 #undef __FUNCT__ 3923 #define __FUNCT__ "TSComputeFunction_Matlab" 3924 /* 3925 TSComputeFunction_Matlab - Calls the function that has been set with 3926 TSSetFunctionMatlab(). 3927 3928 Collective on TS 3929 3930 Input Parameters: 3931 + snes - the TS context 3932 - u - input vector 3933 3934 Output Parameter: 3935 . y - function vector, as set by TSSetFunction() 3936 3937 Notes: 3938 TSComputeFunction() is typically used within nonlinear solvers 3939 implementations, so most users would not generally call this routine 3940 themselves. 3941 3942 Level: developer 3943 3944 .keywords: TS, nonlinear, compute, function 3945 3946 .seealso: TSSetFunction(), TSGetFunction() 3947 */ 3948 PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec u,Vec udot,Vec y, void *ctx) 3949 { 3950 PetscErrorCode ierr; 3951 TSMatlabContext *sctx = (TSMatlabContext *)ctx; 3952 int nlhs = 1,nrhs = 7; 3953 mxArray *plhs[1],*prhs[7]; 3954 long long int lx = 0,lxdot = 0,ly = 0,ls = 0; 3955 3956 PetscFunctionBegin; 3957 PetscValidHeaderSpecific(snes,TS_CLASSID,1); 3958 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 3959 PetscValidHeaderSpecific(udot,VEC_CLASSID,4); 3960 PetscValidHeaderSpecific(y,VEC_CLASSID,5); 3961 PetscCheckSameComm(snes,1,u,3); 3962 PetscCheckSameComm(snes,1,y,5); 3963 3964 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 3965 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 3966 ierr = PetscMemcpy(&lxdot,&udot,sizeof(udot));CHKERRQ(ierr); 3967 ierr = PetscMemcpy(&ly,&y,sizeof(u));CHKERRQ(ierr); 3968 prhs[0] = mxCreateDoubleScalar((double)ls); 3969 prhs[1] = mxCreateDoubleScalar(time); 3970 prhs[2] = mxCreateDoubleScalar((double)lx); 3971 prhs[3] = mxCreateDoubleScalar((double)lxdot); 3972 prhs[4] = mxCreateDoubleScalar((double)ly); 3973 prhs[5] = mxCreateString(sctx->funcname); 3974 prhs[6] = sctx->ctx; 3975 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");CHKERRQ(ierr); 3976 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 3977 mxDestroyArray(prhs[0]); 3978 mxDestroyArray(prhs[1]); 3979 mxDestroyArray(prhs[2]); 3980 mxDestroyArray(prhs[3]); 3981 mxDestroyArray(prhs[4]); 3982 mxDestroyArray(prhs[5]); 3983 mxDestroyArray(plhs[0]); 3984 PetscFunctionReturn(0); 3985 } 3986 3987 3988 #undef __FUNCT__ 3989 #define __FUNCT__ "TSSetFunctionMatlab" 3990 /* 3991 TSSetFunctionMatlab - Sets the function evaluation routine and function 3992 vector for use by the TS routines in solving ODEs 3993 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 3994 3995 Logically Collective on TS 3996 3997 Input Parameters: 3998 + ts - the TS context 3999 - func - function evaluation routine 4000 4001 Calling sequence of func: 4002 $ func (TS ts,PetscReal time,Vec u,Vec udot,Vec f,void *ctx); 4003 4004 Level: beginner 4005 4006 .keywords: TS, nonlinear, set, function 4007 4008 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4009 */ 4010 PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx) 4011 { 4012 PetscErrorCode ierr; 4013 TSMatlabContext *sctx; 4014 4015 PetscFunctionBegin; 4016 /* currently sctx is memory bleed */ 4017 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4018 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4019 /* 4020 This should work, but it doesn't 4021 sctx->ctx = ctx; 4022 mexMakeArrayPersistent(sctx->ctx); 4023 */ 4024 sctx->ctx = mxDuplicateArray(ctx); 4025 ierr = TSSetIFunction(ts,PETSC_NULL,TSComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4026 PetscFunctionReturn(0); 4027 } 4028 4029 #undef __FUNCT__ 4030 #define __FUNCT__ "TSComputeJacobian_Matlab" 4031 /* 4032 TSComputeJacobian_Matlab - Calls the function that has been set with 4033 TSSetJacobianMatlab(). 4034 4035 Collective on TS 4036 4037 Input Parameters: 4038 + ts - the TS context 4039 . u - input vector 4040 . A, B - the matrices 4041 - ctx - user context 4042 4043 Output Parameter: 4044 . flag - structure of the matrix 4045 4046 Level: developer 4047 4048 .keywords: TS, nonlinear, compute, function 4049 4050 .seealso: TSSetFunction(), TSGetFunction() 4051 @*/ 4052 PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec u,Vec udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4053 { 4054 PetscErrorCode ierr; 4055 TSMatlabContext *sctx = (TSMatlabContext *)ctx; 4056 int nlhs = 2,nrhs = 9; 4057 mxArray *plhs[2],*prhs[9]; 4058 long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0; 4059 4060 PetscFunctionBegin; 4061 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4062 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 4063 4064 /* call Matlab function in ctx with arguments u and y */ 4065 4066 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 4067 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4068 ierr = PetscMemcpy(&lxdot,&udot,sizeof(u));CHKERRQ(ierr); 4069 ierr = PetscMemcpy(&lA,A,sizeof(u));CHKERRQ(ierr); 4070 ierr = PetscMemcpy(&lB,B,sizeof(u));CHKERRQ(ierr); 4071 prhs[0] = mxCreateDoubleScalar((double)ls); 4072 prhs[1] = mxCreateDoubleScalar((double)time); 4073 prhs[2] = mxCreateDoubleScalar((double)lx); 4074 prhs[3] = mxCreateDoubleScalar((double)lxdot); 4075 prhs[4] = mxCreateDoubleScalar((double)shift); 4076 prhs[5] = mxCreateDoubleScalar((double)lA); 4077 prhs[6] = mxCreateDoubleScalar((double)lB); 4078 prhs[7] = mxCreateString(sctx->funcname); 4079 prhs[8] = sctx->ctx; 4080 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");CHKERRQ(ierr); 4081 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4082 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4083 mxDestroyArray(prhs[0]); 4084 mxDestroyArray(prhs[1]); 4085 mxDestroyArray(prhs[2]); 4086 mxDestroyArray(prhs[3]); 4087 mxDestroyArray(prhs[4]); 4088 mxDestroyArray(prhs[5]); 4089 mxDestroyArray(prhs[6]); 4090 mxDestroyArray(prhs[7]); 4091 mxDestroyArray(plhs[0]); 4092 mxDestroyArray(plhs[1]); 4093 PetscFunctionReturn(0); 4094 } 4095 4096 4097 #undef __FUNCT__ 4098 #define __FUNCT__ "TSSetJacobianMatlab" 4099 /* 4100 TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4101 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 4102 4103 Logically Collective on TS 4104 4105 Input Parameters: 4106 + ts - the TS context 4107 . A,B - Jacobian matrices 4108 . func - function evaluation routine 4109 - ctx - user context 4110 4111 Calling sequence of func: 4112 $ flag = func (TS ts,PetscReal time,Vec u,Vec udot,Mat A,Mat B,void *ctx); 4113 4114 4115 Level: developer 4116 4117 .keywords: TS, nonlinear, set, function 4118 4119 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4120 */ 4121 PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx) 4122 { 4123 PetscErrorCode ierr; 4124 TSMatlabContext *sctx; 4125 4126 PetscFunctionBegin; 4127 /* currently sctx is memory bleed */ 4128 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4129 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4130 /* 4131 This should work, but it doesn't 4132 sctx->ctx = ctx; 4133 mexMakeArrayPersistent(sctx->ctx); 4134 */ 4135 sctx->ctx = mxDuplicateArray(ctx); 4136 ierr = TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4137 PetscFunctionReturn(0); 4138 } 4139 4140 #undef __FUNCT__ 4141 #define __FUNCT__ "TSMonitor_Matlab" 4142 /* 4143 TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab(). 4144 4145 Collective on TS 4146 4147 .seealso: TSSetFunction(), TSGetFunction() 4148 @*/ 4149 PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec u, void *ctx) 4150 { 4151 PetscErrorCode ierr; 4152 TSMatlabContext *sctx = (TSMatlabContext *)ctx; 4153 int nlhs = 1,nrhs = 6; 4154 mxArray *plhs[1],*prhs[6]; 4155 long long int lx = 0,ls = 0; 4156 4157 PetscFunctionBegin; 4158 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4159 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 4160 4161 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 4162 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4163 prhs[0] = mxCreateDoubleScalar((double)ls); 4164 prhs[1] = mxCreateDoubleScalar((double)it); 4165 prhs[2] = mxCreateDoubleScalar((double)time); 4166 prhs[3] = mxCreateDoubleScalar((double)lx); 4167 prhs[4] = mxCreateString(sctx->funcname); 4168 prhs[5] = sctx->ctx; 4169 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");CHKERRQ(ierr); 4170 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4171 mxDestroyArray(prhs[0]); 4172 mxDestroyArray(prhs[1]); 4173 mxDestroyArray(prhs[2]); 4174 mxDestroyArray(prhs[3]); 4175 mxDestroyArray(prhs[4]); 4176 mxDestroyArray(plhs[0]); 4177 PetscFunctionReturn(0); 4178 } 4179 4180 4181 #undef __FUNCT__ 4182 #define __FUNCT__ "TSMonitorSetMatlab" 4183 /* 4184 TSMonitorSetMatlab - Sets the monitor function from Matlab 4185 4186 Level: developer 4187 4188 .keywords: TS, nonlinear, set, function 4189 4190 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4191 */ 4192 PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx) 4193 { 4194 PetscErrorCode ierr; 4195 TSMatlabContext *sctx; 4196 4197 PetscFunctionBegin; 4198 /* currently sctx is memory bleed */ 4199 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4200 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4201 /* 4202 This should work, but it doesn't 4203 sctx->ctx = ctx; 4204 mexMakeArrayPersistent(sctx->ctx); 4205 */ 4206 sctx->ctx = mxDuplicateArray(ctx); 4207 ierr = TSMonitorSet(ts,TSMonitor_Matlab,sctx,PETSC_NULL);CHKERRQ(ierr); 4208 PetscFunctionReturn(0); 4209 } 4210 #endif 4211 4212 4213 4214 #undef __FUNCT__ 4215 #define __FUNCT__ "TSMonitorLGSolution" 4216 /*@C 4217 TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector 4218 in a time based line graph 4219 4220 Collective on TS 4221 4222 Input Parameters: 4223 + ts - the TS context 4224 . step - current time-step 4225 . ptime - current time 4226 - lg - a line graph object 4227 4228 Level: intermediate 4229 4230 Notes: each process in a parallel run displays its component solutions in a separate window 4231 4232 .keywords: TS, vector, monitor, view 4233 4234 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4235 @*/ 4236 PetscErrorCode TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4237 { 4238 PetscErrorCode ierr; 4239 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 4240 const PetscScalar *yy; 4241 PetscInt dim; 4242 4243 PetscFunctionBegin; 4244 if (!step) { 4245 PetscDrawAxis axis; 4246 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4247 ierr = PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");CHKERRQ(ierr); 4248 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 4249 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 4250 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4251 } 4252 ierr = VecGetArrayRead(u,&yy);CHKERRQ(ierr); 4253 #if defined(PETSC_USE_COMPLEX) 4254 { 4255 PetscReal *yreal; 4256 PetscInt i,n; 4257 ierr = VecGetLocalSize(u,&n);CHKERRQ(ierr); 4258 ierr = PetscMalloc(n*sizeof(PetscReal),&yreal);CHKERRQ(ierr); 4259 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 4260 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 4261 ierr = PetscFree(yreal);CHKERRQ(ierr); 4262 } 4263 #else 4264 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 4265 #endif 4266 ierr = VecRestoreArrayRead(u,&yy);CHKERRQ(ierr); 4267 if (((ctx->howoften > 0) && (!(step % ctx->howoften)) && (step > -1)) || ((ctx->howoften == -1) && (step == -1))){ 4268 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4269 } 4270 PetscFunctionReturn(0); 4271 } 4272 4273 #undef __FUNCT__ 4274 #define __FUNCT__ "TSMonitorLGError" 4275 /*@C 4276 TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the solution vector 4277 in a time based line graph 4278 4279 Collective on TS 4280 4281 Input Parameters: 4282 + ts - the TS context 4283 . step - current time-step 4284 . ptime - current time 4285 - lg - a line graph object 4286 4287 Level: intermediate 4288 4289 Notes: 4290 Only for sequential solves. 4291 4292 The user must provide the solution using TSSetSolutionFunction() to use this monitor. 4293 4294 Options Database Keys: 4295 . -ts_monitor_lg_error - create a graphical monitor of error history 4296 4297 .keywords: TS, vector, monitor, view 4298 4299 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction() 4300 @*/ 4301 PetscErrorCode TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4302 { 4303 PetscErrorCode ierr; 4304 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 4305 const PetscScalar *yy; 4306 Vec y; 4307 PetscInt dim; 4308 4309 PetscFunctionBegin; 4310 if (!step) { 4311 PetscDrawAxis axis; 4312 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4313 ierr = PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Solution");CHKERRQ(ierr); 4314 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 4315 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 4316 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4317 } 4318 ierr = VecDuplicate(u,&y);CHKERRQ(ierr); 4319 ierr = TSComputeSolutionFunction(ts,ptime,y);CHKERRQ(ierr); 4320 ierr = VecAXPY(y,-1.0,u);CHKERRQ(ierr); 4321 ierr = VecGetArrayRead(y,&yy);CHKERRQ(ierr); 4322 #if defined(PETSC_USE_COMPLEX) 4323 { 4324 PetscReal *yreal; 4325 PetscInt i,n; 4326 ierr = VecGetLocalSize(y,&n);CHKERRQ(ierr); 4327 ierr = PetscMalloc(n*sizeof(PetscReal),&yreal);CHKERRQ(ierr); 4328 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 4329 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 4330 ierr = PetscFree(yreal);CHKERRQ(ierr); 4331 } 4332 #else 4333 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 4334 #endif 4335 ierr = VecRestoreArrayRead(y,&yy);CHKERRQ(ierr); 4336 ierr = VecDestroy(&y);CHKERRQ(ierr); 4337 if (((ctx->howoften > 0) && (!(step % ctx->howoften)) && (step > -1)) || ((ctx->howoften == -1) && (step == -1))){ 4338 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4339 } 4340 PetscFunctionReturn(0); 4341 } 4342 4343 #undef __FUNCT__ 4344 #define __FUNCT__ "TSMonitorLGSNESIterations" 4345 PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 4346 { 4347 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 4348 PetscReal x = ptime,y; 4349 PetscErrorCode ierr; 4350 PetscInt its; 4351 4352 PetscFunctionBegin; 4353 if (!n) { 4354 PetscDrawAxis axis; 4355 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4356 ierr = PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");CHKERRQ(ierr); 4357 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4358 ctx->snes_its = 0; 4359 } 4360 ierr = TSGetSNESIterations(ts,&its);CHKERRQ(ierr); 4361 y = its - ctx->snes_its; 4362 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 4363 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))){ 4364 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4365 } 4366 ctx->snes_its = its; 4367 PetscFunctionReturn(0); 4368 } 4369 4370 #undef __FUNCT__ 4371 #define __FUNCT__ "TSMonitorLGKSPIterations" 4372 PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 4373 { 4374 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 4375 PetscReal x = ptime,y; 4376 PetscErrorCode ierr; 4377 PetscInt its; 4378 4379 PetscFunctionBegin; 4380 if (!n) { 4381 PetscDrawAxis axis; 4382 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4383 ierr = PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");CHKERRQ(ierr); 4384 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4385 ctx->ksp_its = 0; 4386 } 4387 ierr = TSGetKSPIterations(ts,&its);CHKERRQ(ierr); 4388 y = its - ctx->ksp_its; 4389 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 4390 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))){ 4391 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4392 } 4393 ctx->ksp_its = its; 4394 PetscFunctionReturn(0); 4395 } 4396 4397 #undef __FUNCT__ 4398 #define __FUNCT__ "TSComputeLinearStability" 4399 /*@ 4400 TSComputeLinearStability - computes the linear stability function at a point 4401 4402 Collective on TS and Vec 4403 4404 Input Parameters: 4405 + ts - the TS context 4406 - xr,xi - real and imaginary part of input arguments 4407 4408 Output Parameters: 4409 . yr,yi - real and imaginary part of function value 4410 4411 Level: developer 4412 4413 .keywords: TS, compute 4414 4415 .seealso: TSSetRHSFunction(), TSComputeIFunction() 4416 @*/ 4417 PetscErrorCode TSComputeLinearStability(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi) 4418 { 4419 PetscErrorCode ierr; 4420 4421 PetscFunctionBegin; 4422 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4423 if (!ts->ops->linearstability) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_SUP,"Linearized stability function not provided for this method"); 4424 ierr = (*ts->ops->linearstability)(ts,xr,xi,yr,yi);CHKERRQ(ierr); 4425 PetscFunctionReturn(0); 4426 } 4427