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 = TSGetAdapt(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 1074 PetscFunctionBegin; 1075 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1076 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 1077 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1078 if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()"); 1079 1080 ierr = PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);CHKERRQ(ierr); 1081 if (classid != TS_FILE_CLASSID) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_ARG_WRONG,"Not TS next in file"); 1082 ierr = PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);CHKERRQ(ierr); 1083 ierr = TSSetType(ts, type);CHKERRQ(ierr); 1084 if (ts->ops->load) { 1085 ierr = (*ts->ops->load)(ts,viewer);CHKERRQ(ierr); 1086 } 1087 ierr = DMCreate(((PetscObject)ts)->comm,&ts->dm);CHKERRQ(ierr); 1088 ierr = DMLoad(ts->dm,viewer);CHKERRQ(ierr); 1089 ierr = DMCreateGlobalVector(ts->dm,&ts->vec_sol);CHKERRQ(ierr); 1090 ierr = VecLoad(ts->vec_sol,viewer);CHKERRQ(ierr); 1091 PetscFunctionReturn(0); 1092 } 1093 1094 #undef __FUNCT__ 1095 #define __FUNCT__ "TSView" 1096 /*@C 1097 TSView - Prints the TS data structure. 1098 1099 Collective on TS 1100 1101 Input Parameters: 1102 + ts - the TS context obtained from TSCreate() 1103 - viewer - visualization context 1104 1105 Options Database Key: 1106 . -ts_view - calls TSView() at end of TSStep() 1107 1108 Notes: 1109 The available visualization contexts include 1110 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 1111 - PETSC_VIEWER_STDOUT_WORLD - synchronized standard 1112 output where only the first processor opens 1113 the file. All other processors send their 1114 data to the first processor to print. 1115 1116 The user can open an alternative visualization context with 1117 PetscViewerASCIIOpen() - output to a specified file. 1118 1119 Level: beginner 1120 1121 .keywords: TS, timestep, view 1122 1123 .seealso: PetscViewerASCIIOpen() 1124 @*/ 1125 PetscErrorCode TSView(TS ts,PetscViewer viewer) 1126 { 1127 PetscErrorCode ierr; 1128 TSType type; 1129 PetscBool iascii,isstring,isundials,isbinary,isdraw; 1130 1131 PetscFunctionBegin; 1132 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1133 if (!viewer) { 1134 ierr = PetscViewerASCIIGetStdout(((PetscObject)ts)->comm,&viewer);CHKERRQ(ierr); 1135 } 1136 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 1137 PetscCheckSameComm(ts,1,viewer,2); 1138 1139 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1140 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 1141 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1142 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1143 if (iascii) { 1144 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)ts,viewer,"TS Object");CHKERRQ(ierr); 1145 ierr = PetscViewerASCIIPrintf(viewer," maximum steps=%D\n",ts->max_steps);CHKERRQ(ierr); 1146 ierr = PetscViewerASCIIPrintf(viewer," maximum time=%G\n",ts->max_time);CHKERRQ(ierr); 1147 if (ts->problem_type == TS_NONLINEAR) { 1148 ierr = PetscViewerASCIIPrintf(viewer," total number of nonlinear solver iterations=%D\n",ts->snes_its);CHKERRQ(ierr); 1149 ierr = PetscViewerASCIIPrintf(viewer," total number of nonlinear solve failures=%D\n",ts->num_snes_failures);CHKERRQ(ierr); 1150 } 1151 ierr = PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",ts->ksp_its);CHKERRQ(ierr); 1152 ierr = PetscViewerASCIIPrintf(viewer," total number of rejected steps=%D\n",ts->reject);CHKERRQ(ierr); 1153 if (ts->ops->view) { 1154 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1155 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1156 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1157 } 1158 } else if (isstring) { 1159 ierr = TSGetType(ts,&type);CHKERRQ(ierr); 1160 ierr = PetscViewerStringSPrintf(viewer," %-7.7s",type);CHKERRQ(ierr); 1161 } else if (isbinary) { 1162 PetscInt classid = TS_FILE_CLASSID; 1163 MPI_Comm comm; 1164 PetscMPIInt rank; 1165 char type[256]; 1166 1167 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 1168 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1169 if (!rank) { 1170 ierr = PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 1171 ierr = PetscStrncpy(type,((PetscObject)ts)->type_name,256);CHKERRQ(ierr); 1172 ierr = PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);CHKERRQ(ierr); 1173 } 1174 if (ts->ops->view) { 1175 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1176 } 1177 ierr = DMView(ts->dm,viewer);CHKERRQ(ierr); 1178 ierr = VecView(ts->vec_sol,viewer);CHKERRQ(ierr); 1179 } else if (isdraw) { 1180 PetscDraw draw; 1181 char str[36]; 1182 PetscReal x,y,bottom,h; 1183 1184 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1185 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 1186 ierr = PetscStrcpy(str,"TS: ");CHKERRQ(ierr); 1187 ierr = PetscStrcat(str,((PetscObject)ts)->type_name);CHKERRQ(ierr); 1188 ierr = PetscDrawBoxedString(draw,x,y,PETSC_DRAW_BLACK,PETSC_DRAW_BLACK,str,PETSC_NULL,&h);CHKERRQ(ierr); 1189 bottom = y - h; 1190 ierr = PetscDrawPushCurrentPoint(draw,x,bottom);CHKERRQ(ierr); 1191 if (ts->ops->view) { 1192 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1193 } 1194 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 1195 } 1196 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1197 ierr = PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&isundials);CHKERRQ(ierr); 1198 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1199 PetscFunctionReturn(0); 1200 } 1201 1202 1203 #undef __FUNCT__ 1204 #define __FUNCT__ "TSSetApplicationContext" 1205 /*@ 1206 TSSetApplicationContext - Sets an optional user-defined context for 1207 the timesteppers. 1208 1209 Logically Collective on TS 1210 1211 Input Parameters: 1212 + ts - the TS context obtained from TSCreate() 1213 - usrP - optional user context 1214 1215 Level: intermediate 1216 1217 .keywords: TS, timestep, set, application, context 1218 1219 .seealso: TSGetApplicationContext() 1220 @*/ 1221 PetscErrorCode TSSetApplicationContext(TS ts,void *usrP) 1222 { 1223 PetscFunctionBegin; 1224 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1225 ts->user = usrP; 1226 PetscFunctionReturn(0); 1227 } 1228 1229 #undef __FUNCT__ 1230 #define __FUNCT__ "TSGetApplicationContext" 1231 /*@ 1232 TSGetApplicationContext - Gets the user-defined context for the 1233 timestepper. 1234 1235 Not Collective 1236 1237 Input Parameter: 1238 . ts - the TS context obtained from TSCreate() 1239 1240 Output Parameter: 1241 . usrP - user context 1242 1243 Level: intermediate 1244 1245 .keywords: TS, timestep, get, application, context 1246 1247 .seealso: TSSetApplicationContext() 1248 @*/ 1249 PetscErrorCode TSGetApplicationContext(TS ts,void *usrP) 1250 { 1251 PetscFunctionBegin; 1252 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1253 *(void**)usrP = ts->user; 1254 PetscFunctionReturn(0); 1255 } 1256 1257 #undef __FUNCT__ 1258 #define __FUNCT__ "TSGetTimeStepNumber" 1259 /*@ 1260 TSGetTimeStepNumber - Gets the number of time steps completed. 1261 1262 Not Collective 1263 1264 Input Parameter: 1265 . ts - the TS context obtained from TSCreate() 1266 1267 Output Parameter: 1268 . iter - number of steps completed so far 1269 1270 Level: intermediate 1271 1272 .keywords: TS, timestep, get, iteration, number 1273 .seealso: TSGetTime(), TSGetTimeStep(), TSSetPreStep(), TSSetPreStage(), TSSetPostStep() 1274 @*/ 1275 PetscErrorCode TSGetTimeStepNumber(TS ts,PetscInt* iter) 1276 { 1277 PetscFunctionBegin; 1278 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1279 PetscValidIntPointer(iter,2); 1280 *iter = ts->steps; 1281 PetscFunctionReturn(0); 1282 } 1283 1284 #undef __FUNCT__ 1285 #define __FUNCT__ "TSSetInitialTimeStep" 1286 /*@ 1287 TSSetInitialTimeStep - Sets the initial timestep to be used, 1288 as well as the initial time. 1289 1290 Logically Collective on TS 1291 1292 Input Parameters: 1293 + ts - the TS context obtained from TSCreate() 1294 . initial_time - the initial time 1295 - time_step - the size of the timestep 1296 1297 Level: intermediate 1298 1299 .seealso: TSSetTimeStep(), TSGetTimeStep() 1300 1301 .keywords: TS, set, initial, timestep 1302 @*/ 1303 PetscErrorCode TSSetInitialTimeStep(TS ts,PetscReal initial_time,PetscReal time_step) 1304 { 1305 PetscErrorCode ierr; 1306 1307 PetscFunctionBegin; 1308 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1309 ierr = TSSetTimeStep(ts,time_step);CHKERRQ(ierr); 1310 ierr = TSSetTime(ts,initial_time);CHKERRQ(ierr); 1311 PetscFunctionReturn(0); 1312 } 1313 1314 #undef __FUNCT__ 1315 #define __FUNCT__ "TSSetTimeStep" 1316 /*@ 1317 TSSetTimeStep - Allows one to reset the timestep at any time, 1318 useful for simple pseudo-timestepping codes. 1319 1320 Logically Collective on TS 1321 1322 Input Parameters: 1323 + ts - the TS context obtained from TSCreate() 1324 - time_step - the size of the timestep 1325 1326 Level: intermediate 1327 1328 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 1329 1330 .keywords: TS, set, timestep 1331 @*/ 1332 PetscErrorCode TSSetTimeStep(TS ts,PetscReal time_step) 1333 { 1334 PetscFunctionBegin; 1335 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1336 PetscValidLogicalCollectiveReal(ts,time_step,2); 1337 ts->time_step = time_step; 1338 ts->time_step_orig = time_step; 1339 PetscFunctionReturn(0); 1340 } 1341 1342 #undef __FUNCT__ 1343 #define __FUNCT__ "TSSetExactFinalTime" 1344 /*@ 1345 TSSetExactFinalTime - Determines whether to adapt the final time step to 1346 match the exact final time, interpolate solution to the exact final time, 1347 or just return at the final time TS computed. 1348 1349 Logically Collective on TS 1350 1351 Input Parameter: 1352 + ts - the time-step context 1353 - eftopt - exact final time option 1354 1355 Level: beginner 1356 1357 .seealso: TSExactFinalTime 1358 @*/ 1359 PetscErrorCode TSSetExactFinalTime(TS ts,TSExactFinalTimeOption eftopt) 1360 { 1361 PetscFunctionBegin; 1362 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1363 PetscValidLogicalCollectiveEnum(ts,eftopt,2); 1364 ts->exact_final_time = eftopt; 1365 PetscFunctionReturn(0); 1366 } 1367 1368 #undef __FUNCT__ 1369 #define __FUNCT__ "TSGetTimeStep" 1370 /*@ 1371 TSGetTimeStep - Gets the current timestep size. 1372 1373 Not Collective 1374 1375 Input Parameter: 1376 . ts - the TS context obtained from TSCreate() 1377 1378 Output Parameter: 1379 . dt - the current timestep size 1380 1381 Level: intermediate 1382 1383 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 1384 1385 .keywords: TS, get, timestep 1386 @*/ 1387 PetscErrorCode TSGetTimeStep(TS ts,PetscReal* dt) 1388 { 1389 PetscFunctionBegin; 1390 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1391 PetscValidRealPointer(dt,2); 1392 *dt = ts->time_step; 1393 PetscFunctionReturn(0); 1394 } 1395 1396 #undef __FUNCT__ 1397 #define __FUNCT__ "TSGetSolution" 1398 /*@ 1399 TSGetSolution - Returns the solution at the present timestep. It 1400 is valid to call this routine inside the function that you are evaluating 1401 in order to move to the new timestep. This vector not changed until 1402 the solution at the next timestep has been calculated. 1403 1404 Not Collective, but Vec returned is parallel if TS is parallel 1405 1406 Input Parameter: 1407 . ts - the TS context obtained from TSCreate() 1408 1409 Output Parameter: 1410 . v - the vector containing the solution 1411 1412 Level: intermediate 1413 1414 .seealso: TSGetTimeStep() 1415 1416 .keywords: TS, timestep, get, solution 1417 @*/ 1418 PetscErrorCode TSGetSolution(TS ts,Vec *v) 1419 { 1420 PetscFunctionBegin; 1421 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1422 PetscValidPointer(v,2); 1423 *v = ts->vec_sol; 1424 PetscFunctionReturn(0); 1425 } 1426 1427 /* ----- Routines to initialize and destroy a timestepper ---- */ 1428 #undef __FUNCT__ 1429 #define __FUNCT__ "TSSetProblemType" 1430 /*@ 1431 TSSetProblemType - Sets the type of problem to be solved. 1432 1433 Not collective 1434 1435 Input Parameters: 1436 + ts - The TS 1437 - type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms 1438 .vb 1439 U_t - A U = 0 (linear) 1440 U_t - A(t) U = 0 (linear) 1441 F(t,U,U_t) = 0 (nonlinear) 1442 .ve 1443 1444 Level: beginner 1445 1446 .keywords: TS, problem type 1447 .seealso: TSSetUp(), TSProblemType, TS 1448 @*/ 1449 PetscErrorCode TSSetProblemType(TS ts, TSProblemType type) 1450 { 1451 PetscErrorCode ierr; 1452 1453 PetscFunctionBegin; 1454 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1455 ts->problem_type = type; 1456 if (type == TS_LINEAR) { 1457 SNES snes; 1458 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1459 ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 1460 } 1461 PetscFunctionReturn(0); 1462 } 1463 1464 #undef __FUNCT__ 1465 #define __FUNCT__ "TSGetProblemType" 1466 /*@C 1467 TSGetProblemType - Gets the type of problem to be solved. 1468 1469 Not collective 1470 1471 Input Parameter: 1472 . ts - The TS 1473 1474 Output Parameter: 1475 . type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms 1476 .vb 1477 M U_t = A U 1478 M(t) U_t = A(t) U 1479 F(t,U,U_t) 1480 .ve 1481 1482 Level: beginner 1483 1484 .keywords: TS, problem type 1485 .seealso: TSSetUp(), TSProblemType, TS 1486 @*/ 1487 PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type) 1488 { 1489 PetscFunctionBegin; 1490 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1491 PetscValidIntPointer(type,2); 1492 *type = ts->problem_type; 1493 PetscFunctionReturn(0); 1494 } 1495 1496 #undef __FUNCT__ 1497 #define __FUNCT__ "TSSetUp" 1498 /*@ 1499 TSSetUp - Sets up the internal data structures for the later use 1500 of a timestepper. 1501 1502 Collective on TS 1503 1504 Input Parameter: 1505 . ts - the TS context obtained from TSCreate() 1506 1507 Notes: 1508 For basic use of the TS solvers the user need not explicitly call 1509 TSSetUp(), since these actions will automatically occur during 1510 the call to TSStep(). However, if one wishes to control this 1511 phase separately, TSSetUp() should be called after TSCreate() 1512 and optional routines of the form TSSetXXX(), but before TSStep(). 1513 1514 Level: advanced 1515 1516 .keywords: TS, timestep, setup 1517 1518 .seealso: TSCreate(), TSStep(), TSDestroy() 1519 @*/ 1520 PetscErrorCode TSSetUp(TS ts) 1521 { 1522 PetscErrorCode ierr; 1523 DM dm; 1524 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 1525 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 1526 TSIJacobian ijac; 1527 TSRHSJacobian rhsjac; 1528 1529 PetscFunctionBegin; 1530 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1531 if (ts->setupcalled) PetscFunctionReturn(0); 1532 1533 if (!((PetscObject)ts)->type_name) { 1534 ierr = TSSetType(ts,TSEULER);CHKERRQ(ierr); 1535 } 1536 1537 if (!ts->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetSolution() first"); 1538 1539 ierr = TSGetAdapt(ts,&ts->adapt);CHKERRQ(ierr); 1540 1541 if (ts->ops->setup) { 1542 ierr = (*ts->ops->setup)(ts);CHKERRQ(ierr); 1543 } 1544 1545 /* in the case where we've set a DMTSFunction or what have you, we need the default SNESFunction 1546 to be set right but can't do it elsewhere due to the overreliance on ctx=ts. 1547 */ 1548 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1549 ierr = DMSNESGetFunction(dm,&func,PETSC_NULL);CHKERRQ(ierr); 1550 if (!func) { 1551 ierr =DMSNESSetFunction(dm,SNESTSFormFunction,ts);CHKERRQ(ierr); 1552 } 1553 /* if the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it. 1554 Otherwise, the SNES will use coloring internally to form the Jacobian. 1555 */ 1556 ierr = DMSNESGetJacobian(dm,&jac,PETSC_NULL);CHKERRQ(ierr); 1557 ierr = DMTSGetIJacobian(dm,&ijac,PETSC_NULL);CHKERRQ(ierr); 1558 ierr = DMTSGetRHSJacobian(dm,&rhsjac,PETSC_NULL);CHKERRQ(ierr); 1559 if (!jac && (ijac || rhsjac)) { 1560 ierr = DMSNESSetJacobian(dm,SNESTSFormJacobian,ts);CHKERRQ(ierr); 1561 } 1562 ts->setupcalled = PETSC_TRUE; 1563 PetscFunctionReturn(0); 1564 } 1565 1566 #undef __FUNCT__ 1567 #define __FUNCT__ "TSReset" 1568 /*@ 1569 TSReset - Resets a TS context and removes any allocated Vecs and Mats. 1570 1571 Collective on TS 1572 1573 Input Parameter: 1574 . ts - the TS context obtained from TSCreate() 1575 1576 Level: beginner 1577 1578 .keywords: TS, timestep, reset 1579 1580 .seealso: TSCreate(), TSSetup(), TSDestroy() 1581 @*/ 1582 PetscErrorCode TSReset(TS ts) 1583 { 1584 PetscErrorCode ierr; 1585 1586 PetscFunctionBegin; 1587 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1588 if (ts->ops->reset) { 1589 ierr = (*ts->ops->reset)(ts);CHKERRQ(ierr); 1590 } 1591 if (ts->snes) {ierr = SNESReset(ts->snes);CHKERRQ(ierr);} 1592 ierr = MatDestroy(&ts->Arhs);CHKERRQ(ierr); 1593 ierr = MatDestroy(&ts->Brhs);CHKERRQ(ierr); 1594 ierr = VecDestroy(&ts->Frhs);CHKERRQ(ierr); 1595 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 1596 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 1597 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 1598 ierr = VecDestroyVecs(ts->nwork,&ts->work);CHKERRQ(ierr); 1599 ts->setupcalled = PETSC_FALSE; 1600 PetscFunctionReturn(0); 1601 } 1602 1603 #undef __FUNCT__ 1604 #define __FUNCT__ "TSDestroy" 1605 /*@ 1606 TSDestroy - Destroys the timestepper context that was created 1607 with TSCreate(). 1608 1609 Collective on TS 1610 1611 Input Parameter: 1612 . ts - the TS context obtained from TSCreate() 1613 1614 Level: beginner 1615 1616 .keywords: TS, timestepper, destroy 1617 1618 .seealso: TSCreate(), TSSetUp(), TSSolve() 1619 @*/ 1620 PetscErrorCode TSDestroy(TS *ts) 1621 { 1622 PetscErrorCode ierr; 1623 1624 PetscFunctionBegin; 1625 if (!*ts) PetscFunctionReturn(0); 1626 PetscValidHeaderSpecific((*ts),TS_CLASSID,1); 1627 if (--((PetscObject)(*ts))->refct > 0) {*ts = 0; PetscFunctionReturn(0);} 1628 1629 ierr = TSReset((*ts));CHKERRQ(ierr); 1630 1631 /* if memory was published with AMS then destroy it */ 1632 ierr = PetscObjectDepublish((*ts));CHKERRQ(ierr); 1633 if ((*ts)->ops->destroy) {ierr = (*(*ts)->ops->destroy)((*ts));CHKERRQ(ierr);} 1634 1635 ierr = TSAdaptDestroy(&(*ts)->adapt);CHKERRQ(ierr); 1636 ierr = SNESDestroy(&(*ts)->snes);CHKERRQ(ierr); 1637 ierr = DMDestroy(&(*ts)->dm);CHKERRQ(ierr); 1638 ierr = TSMonitorCancel((*ts));CHKERRQ(ierr); 1639 1640 ierr = PetscHeaderDestroy(ts);CHKERRQ(ierr); 1641 PetscFunctionReturn(0); 1642 } 1643 1644 #undef __FUNCT__ 1645 #define __FUNCT__ "TSGetSNES" 1646 /*@ 1647 TSGetSNES - Returns the SNES (nonlinear solver) associated with 1648 a TS (timestepper) context. Valid only for nonlinear problems. 1649 1650 Not Collective, but SNES is parallel if TS is parallel 1651 1652 Input Parameter: 1653 . ts - the TS context obtained from TSCreate() 1654 1655 Output Parameter: 1656 . snes - the nonlinear solver context 1657 1658 Notes: 1659 The user can then directly manipulate the SNES context to set various 1660 options, etc. Likewise, the user can then extract and manipulate the 1661 KSP, KSP, and PC contexts as well. 1662 1663 TSGetSNES() does not work for integrators that do not use SNES; in 1664 this case TSGetSNES() returns PETSC_NULL in snes. 1665 1666 Level: beginner 1667 1668 .keywords: timestep, get, SNES 1669 @*/ 1670 PetscErrorCode TSGetSNES(TS ts,SNES *snes) 1671 { 1672 PetscErrorCode ierr; 1673 1674 PetscFunctionBegin; 1675 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1676 PetscValidPointer(snes,2); 1677 if (!ts->snes) { 1678 ierr = SNESCreate(((PetscObject)ts)->comm,&ts->snes);CHKERRQ(ierr); 1679 ierr = SNESSetFunction(ts->snes,PETSC_NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 1680 ierr = PetscLogObjectParent(ts,ts->snes);CHKERRQ(ierr); 1681 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);CHKERRQ(ierr); 1682 if (ts->dm) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 1683 if (ts->problem_type == TS_LINEAR) { 1684 ierr = SNESSetType(ts->snes,SNESKSPONLY);CHKERRQ(ierr); 1685 } 1686 } 1687 *snes = ts->snes; 1688 PetscFunctionReturn(0); 1689 } 1690 1691 #undef __FUNCT__ 1692 #define __FUNCT__ "TSGetKSP" 1693 /*@ 1694 TSGetKSP - Returns the KSP (linear solver) associated with 1695 a TS (timestepper) context. 1696 1697 Not Collective, but KSP is parallel if TS is parallel 1698 1699 Input Parameter: 1700 . ts - the TS context obtained from TSCreate() 1701 1702 Output Parameter: 1703 . ksp - the nonlinear solver context 1704 1705 Notes: 1706 The user can then directly manipulate the KSP context to set various 1707 options, etc. Likewise, the user can then extract and manipulate the 1708 KSP and PC contexts as well. 1709 1710 TSGetKSP() does not work for integrators that do not use KSP; 1711 in this case TSGetKSP() returns PETSC_NULL in ksp. 1712 1713 Level: beginner 1714 1715 .keywords: timestep, get, KSP 1716 @*/ 1717 PetscErrorCode TSGetKSP(TS ts,KSP *ksp) 1718 { 1719 PetscErrorCode ierr; 1720 SNES snes; 1721 1722 PetscFunctionBegin; 1723 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1724 PetscValidPointer(ksp,2); 1725 if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first"); 1726 if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()"); 1727 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1728 ierr = SNESGetKSP(snes,ksp);CHKERRQ(ierr); 1729 PetscFunctionReturn(0); 1730 } 1731 1732 /* ----------- Routines to set solver parameters ---------- */ 1733 1734 #undef __FUNCT__ 1735 #define __FUNCT__ "TSGetDuration" 1736 /*@ 1737 TSGetDuration - Gets the maximum number of timesteps to use and 1738 maximum time for iteration. 1739 1740 Not Collective 1741 1742 Input Parameters: 1743 + ts - the TS context obtained from TSCreate() 1744 . maxsteps - maximum number of iterations to use, or PETSC_NULL 1745 - maxtime - final time to iterate to, or PETSC_NULL 1746 1747 Level: intermediate 1748 1749 .keywords: TS, timestep, get, maximum, iterations, time 1750 @*/ 1751 PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime) 1752 { 1753 PetscFunctionBegin; 1754 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1755 if (maxsteps) { 1756 PetscValidIntPointer(maxsteps,2); 1757 *maxsteps = ts->max_steps; 1758 } 1759 if (maxtime) { 1760 PetscValidScalarPointer(maxtime,3); 1761 *maxtime = ts->max_time; 1762 } 1763 PetscFunctionReturn(0); 1764 } 1765 1766 #undef __FUNCT__ 1767 #define __FUNCT__ "TSSetDuration" 1768 /*@ 1769 TSSetDuration - Sets the maximum number of timesteps to use and 1770 maximum time for iteration. 1771 1772 Logically Collective on TS 1773 1774 Input Parameters: 1775 + ts - the TS context obtained from TSCreate() 1776 . maxsteps - maximum number of iterations to use 1777 - maxtime - final time to iterate to 1778 1779 Options Database Keys: 1780 . -ts_max_steps <maxsteps> - Sets maxsteps 1781 . -ts_final_time <maxtime> - Sets maxtime 1782 1783 Notes: 1784 The default maximum number of iterations is 5000. Default time is 5.0 1785 1786 Level: intermediate 1787 1788 .keywords: TS, timestep, set, maximum, iterations 1789 1790 .seealso: TSSetExactFinalTime() 1791 @*/ 1792 PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime) 1793 { 1794 PetscFunctionBegin; 1795 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1796 PetscValidLogicalCollectiveInt(ts,maxsteps,2); 1797 PetscValidLogicalCollectiveReal(ts,maxtime,2); 1798 if (maxsteps >= 0) ts->max_steps = maxsteps; 1799 if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime; 1800 PetscFunctionReturn(0); 1801 } 1802 1803 #undef __FUNCT__ 1804 #define __FUNCT__ "TSSetSolution" 1805 /*@ 1806 TSSetSolution - Sets the initial solution vector 1807 for use by the TS routines. 1808 1809 Logically Collective on TS and Vec 1810 1811 Input Parameters: 1812 + ts - the TS context obtained from TSCreate() 1813 - u - the solution vector 1814 1815 Level: beginner 1816 1817 .keywords: TS, timestep, set, solution, initial conditions 1818 @*/ 1819 PetscErrorCode TSSetSolution(TS ts,Vec u) 1820 { 1821 PetscErrorCode ierr; 1822 DM dm; 1823 1824 PetscFunctionBegin; 1825 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1826 PetscValidHeaderSpecific(u,VEC_CLASSID,2); 1827 ierr = PetscObjectReference((PetscObject)u);CHKERRQ(ierr); 1828 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 1829 ts->vec_sol = u; 1830 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1831 ierr = DMShellSetGlobalVector(dm,u);CHKERRQ(ierr); 1832 PetscFunctionReturn(0); 1833 } 1834 1835 #undef __FUNCT__ 1836 #define __FUNCT__ "TSSetPreStep" 1837 /*@C 1838 TSSetPreStep - Sets the general-purpose function 1839 called once at the beginning of each time step. 1840 1841 Logically Collective on TS 1842 1843 Input Parameters: 1844 + ts - The TS context obtained from TSCreate() 1845 - func - The function 1846 1847 Calling sequence of func: 1848 . func (TS ts); 1849 1850 Level: intermediate 1851 1852 Note: 1853 If a step is rejected, TSStep() will call this routine again before each attempt. 1854 The last completed time step number can be queried using TSGetTimeStepNumber(), the 1855 size of the step being attempted can be obtained using TSGetTimeStep(). 1856 1857 .keywords: TS, timestep 1858 .seealso: TSSetPreStage(), TSSetPostStep(), TSStep() 1859 @*/ 1860 PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS)) 1861 { 1862 PetscFunctionBegin; 1863 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1864 ts->ops->prestep = func; 1865 PetscFunctionReturn(0); 1866 } 1867 1868 #undef __FUNCT__ 1869 #define __FUNCT__ "TSPreStep" 1870 /*@ 1871 TSPreStep - Runs the user-defined pre-step function. 1872 1873 Collective on TS 1874 1875 Input Parameters: 1876 . ts - The TS context obtained from TSCreate() 1877 1878 Notes: 1879 TSPreStep() is typically used within time stepping implementations, 1880 so most users would not generally call this routine themselves. 1881 1882 Level: developer 1883 1884 .keywords: TS, timestep 1885 .seealso: TSSetPreStep(), TSPreStage(), TSPostStep() 1886 @*/ 1887 PetscErrorCode TSPreStep(TS ts) 1888 { 1889 PetscErrorCode ierr; 1890 1891 PetscFunctionBegin; 1892 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1893 if (ts->ops->prestep) { 1894 PetscStackPush("TS PreStep function"); 1895 ierr = (*ts->ops->prestep)(ts);CHKERRQ(ierr); 1896 PetscStackPop; 1897 } 1898 PetscFunctionReturn(0); 1899 } 1900 1901 #undef __FUNCT__ 1902 #define __FUNCT__ "TSSetPreStage" 1903 /*@C 1904 TSSetPreStage - Sets the general-purpose function 1905 called once at the beginning of each stage. 1906 1907 Logically Collective on TS 1908 1909 Input Parameters: 1910 + ts - The TS context obtained from TSCreate() 1911 - func - The function 1912 1913 Calling sequence of func: 1914 . PetscErrorCode func(TS ts, PetscReal stagetime); 1915 1916 Level: intermediate 1917 1918 Note: 1919 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 1920 The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being 1921 attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime(). 1922 1923 .keywords: TS, timestep 1924 .seealso: TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 1925 @*/ 1926 PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS,PetscReal)) 1927 { 1928 PetscFunctionBegin; 1929 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1930 ts->ops->prestage = func; 1931 PetscFunctionReturn(0); 1932 } 1933 1934 #undef __FUNCT__ 1935 #define __FUNCT__ "TSPreStage" 1936 /*@ 1937 TSPreStage - Runs the user-defined pre-stage function set using TSSetPreStage() 1938 1939 Collective on TS 1940 1941 Input Parameters: 1942 . ts - The TS context obtained from TSCreate() 1943 1944 Notes: 1945 TSPreStage() is typically used within time stepping implementations, 1946 most users would not generally call this routine themselves. 1947 1948 Level: developer 1949 1950 .keywords: TS, timestep 1951 .seealso: TSSetPreStep(), TSPreStep(), TSPostStep() 1952 @*/ 1953 PetscErrorCode TSPreStage(TS ts, PetscReal stagetime) 1954 { 1955 PetscErrorCode ierr; 1956 1957 PetscFunctionBegin; 1958 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1959 if (ts->ops->prestage) { 1960 PetscStackPush("TS PreStage function"); 1961 ierr = (*ts->ops->prestage)(ts,stagetime);CHKERRQ(ierr); 1962 PetscStackPop; 1963 } 1964 PetscFunctionReturn(0); 1965 } 1966 1967 #undef __FUNCT__ 1968 #define __FUNCT__ "TSSetPostStep" 1969 /*@C 1970 TSSetPostStep - Sets the general-purpose function 1971 called once at the end of each time step. 1972 1973 Logically Collective on TS 1974 1975 Input Parameters: 1976 + ts - The TS context obtained from TSCreate() 1977 - func - The function 1978 1979 Calling sequence of func: 1980 $ func (TS ts); 1981 1982 Level: intermediate 1983 1984 .keywords: TS, timestep 1985 .seealso: TSSetPreStep(), TSSetPreStage(), TSGetTimeStep(), TSGetTimeStepNumber(), TSGetTime() 1986 @*/ 1987 PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS)) 1988 { 1989 PetscFunctionBegin; 1990 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1991 ts->ops->poststep = func; 1992 PetscFunctionReturn(0); 1993 } 1994 1995 #undef __FUNCT__ 1996 #define __FUNCT__ "TSPostStep" 1997 /*@ 1998 TSPostStep - Runs the user-defined post-step function. 1999 2000 Collective on TS 2001 2002 Input Parameters: 2003 . ts - The TS context obtained from TSCreate() 2004 2005 Notes: 2006 TSPostStep() is typically used within time stepping implementations, 2007 so most users would not generally call this routine themselves. 2008 2009 Level: developer 2010 2011 .keywords: TS, timestep 2012 @*/ 2013 PetscErrorCode TSPostStep(TS ts) 2014 { 2015 PetscErrorCode ierr; 2016 2017 PetscFunctionBegin; 2018 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2019 if (ts->ops->poststep) { 2020 PetscStackPush("TS PostStep function"); 2021 ierr = (*ts->ops->poststep)(ts);CHKERRQ(ierr); 2022 PetscStackPop; 2023 } 2024 PetscFunctionReturn(0); 2025 } 2026 2027 /* ------------ Routines to set performance monitoring options ----------- */ 2028 2029 #undef __FUNCT__ 2030 #define __FUNCT__ "TSMonitorSet" 2031 /*@C 2032 TSMonitorSet - Sets an ADDITIONAL function that is to be used at every 2033 timestep to display the iteration's progress. 2034 2035 Logically Collective on TS 2036 2037 Input Parameters: 2038 + ts - the TS context obtained from TSCreate() 2039 . monitor - monitoring routine 2040 . mctx - [optional] user-defined context for private data for the 2041 monitor routine (use PETSC_NULL if no context is desired) 2042 - monitordestroy - [optional] routine that frees monitor context 2043 (may be PETSC_NULL) 2044 2045 Calling sequence of monitor: 2046 $ int monitor(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx) 2047 2048 + ts - the TS context 2049 . 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 2050 been interpolated to) 2051 . time - current time 2052 . u - current iterate 2053 - mctx - [optional] monitoring context 2054 2055 Notes: 2056 This routine adds an additional monitor to the list of monitors that 2057 already has been loaded. 2058 2059 Fortran notes: Only a single monitor function can be set for each TS object 2060 2061 Level: intermediate 2062 2063 .keywords: TS, timestep, set, monitor 2064 2065 .seealso: TSMonitorDefault(), TSMonitorCancel() 2066 @*/ 2067 PetscErrorCode TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**)) 2068 { 2069 PetscFunctionBegin; 2070 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2071 if (ts->numbermonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 2072 ts->monitor[ts->numbermonitors] = monitor; 2073 ts->monitordestroy[ts->numbermonitors] = mdestroy; 2074 ts->monitorcontext[ts->numbermonitors++] = (void*)mctx; 2075 PetscFunctionReturn(0); 2076 } 2077 2078 #undef __FUNCT__ 2079 #define __FUNCT__ "TSMonitorCancel" 2080 /*@C 2081 TSMonitorCancel - Clears all the monitors that have been set on a time-step object. 2082 2083 Logically Collective on TS 2084 2085 Input Parameters: 2086 . ts - the TS context obtained from TSCreate() 2087 2088 Notes: 2089 There is no way to remove a single, specific monitor. 2090 2091 Level: intermediate 2092 2093 .keywords: TS, timestep, set, monitor 2094 2095 .seealso: TSMonitorDefault(), TSMonitorSet() 2096 @*/ 2097 PetscErrorCode TSMonitorCancel(TS ts) 2098 { 2099 PetscErrorCode ierr; 2100 PetscInt i; 2101 2102 PetscFunctionBegin; 2103 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2104 for (i=0; i<ts->numbermonitors; i++) { 2105 if (ts->monitordestroy[i]) { 2106 ierr = (*ts->monitordestroy[i])(&ts->monitorcontext[i]);CHKERRQ(ierr); 2107 } 2108 } 2109 ts->numbermonitors = 0; 2110 PetscFunctionReturn(0); 2111 } 2112 2113 #undef __FUNCT__ 2114 #define __FUNCT__ "TSMonitorDefault" 2115 /*@ 2116 TSMonitorDefault - Sets the Default monitor 2117 2118 Level: intermediate 2119 2120 .keywords: TS, set, monitor 2121 2122 .seealso: TSMonitorDefault(), TSMonitorSet() 2123 @*/ 2124 PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy) 2125 { 2126 PetscErrorCode ierr; 2127 PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(((PetscObject)ts)->comm); 2128 2129 PetscFunctionBegin; 2130 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 2131 ierr = PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g\n",step,(double)ts->time_step,(double)ptime);CHKERRQ(ierr); 2132 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 2133 PetscFunctionReturn(0); 2134 } 2135 2136 #undef __FUNCT__ 2137 #define __FUNCT__ "TSSetRetainStages" 2138 /*@ 2139 TSSetRetainStages - Request that all stages in the upcoming step be stored so that interpolation will be available. 2140 2141 Logically Collective on TS 2142 2143 Input Argument: 2144 . ts - time stepping context 2145 2146 Output Argument: 2147 . flg - PETSC_TRUE or PETSC_FALSE 2148 2149 Level: intermediate 2150 2151 .keywords: TS, set 2152 2153 .seealso: TSInterpolate(), TSSetPostStep() 2154 @*/ 2155 PetscErrorCode TSSetRetainStages(TS ts,PetscBool flg) 2156 { 2157 PetscFunctionBegin; 2158 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2159 ts->retain_stages = flg; 2160 PetscFunctionReturn(0); 2161 } 2162 2163 #undef __FUNCT__ 2164 #define __FUNCT__ "TSInterpolate" 2165 /*@ 2166 TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval 2167 2168 Collective on TS 2169 2170 Input Argument: 2171 + ts - time stepping context 2172 - t - time to interpolate to 2173 2174 Output Argument: 2175 . U - state at given time 2176 2177 Notes: 2178 The user should call TSSetRetainStages() before taking a step in which interpolation will be requested. 2179 2180 Level: intermediate 2181 2182 Developer Notes: 2183 TSInterpolate() and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints. 2184 2185 .keywords: TS, set 2186 2187 .seealso: TSSetRetainStages(), TSSetPostStep() 2188 @*/ 2189 PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec U) 2190 { 2191 PetscErrorCode ierr; 2192 2193 PetscFunctionBegin; 2194 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2195 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); 2196 if (!ts->ops->interpolate) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name); 2197 ierr = (*ts->ops->interpolate)(ts,t,U);CHKERRQ(ierr); 2198 PetscFunctionReturn(0); 2199 } 2200 2201 #undef __FUNCT__ 2202 #define __FUNCT__ "TSStep" 2203 /*@ 2204 TSStep - Steps one time step 2205 2206 Collective on TS 2207 2208 Input Parameter: 2209 . ts - the TS context obtained from TSCreate() 2210 2211 Level: intermediate 2212 2213 Notes: 2214 The hook set using TSSetPreStep() is called before each attempt to take the step. In general, the time step size may 2215 be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages. 2216 2217 This may over-step the final time provided in TSSetDuration() depending on the time-step used. TSSolve() interpolates to exactly the 2218 time provided in TSSetDuration(). One can use TSInterpolate() to determine an interpolated solution within the final timestep. 2219 2220 .keywords: TS, timestep, solve 2221 2222 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSInterpolate() 2223 @*/ 2224 PetscErrorCode TSStep(TS ts) 2225 { 2226 PetscReal ptime_prev; 2227 PetscErrorCode ierr; 2228 2229 PetscFunctionBegin; 2230 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2231 ierr = TSSetUp(ts);CHKERRQ(ierr); 2232 2233 ts->reason = TS_CONVERGED_ITERATING; 2234 2235 ptime_prev = ts->ptime; 2236 ierr = PetscLogEventBegin(TS_Step,ts,0,0,0);CHKERRQ(ierr); 2237 ierr = (*ts->ops->step)(ts);CHKERRQ(ierr); 2238 ierr = PetscLogEventEnd(TS_Step,ts,0,0,0);CHKERRQ(ierr); 2239 ts->time_step_prev = ts->ptime - ptime_prev; 2240 2241 if (ts->reason < 0) { 2242 if (ts->errorifstepfailed) { 2243 if (ts->reason == TS_DIVERGED_NONLINEAR_SOLVE) { 2244 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]); 2245 } else SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]); 2246 } 2247 } else if (!ts->reason) { 2248 if (ts->steps >= ts->max_steps) 2249 ts->reason = TS_CONVERGED_ITS; 2250 else if (ts->ptime >= ts->max_time) 2251 ts->reason = TS_CONVERGED_TIME; 2252 } 2253 2254 PetscFunctionReturn(0); 2255 } 2256 2257 #undef __FUNCT__ 2258 #define __FUNCT__ "TSEvaluateStep" 2259 /*@ 2260 TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy. 2261 2262 Collective on TS 2263 2264 Input Arguments: 2265 + ts - time stepping context 2266 . order - desired order of accuracy 2267 - done - whether the step was evaluated at this order (pass PETSC_NULL to generate an error if not available) 2268 2269 Output Arguments: 2270 . U - state at the end of the current step 2271 2272 Level: advanced 2273 2274 Notes: 2275 This function cannot be called until all stages have been evaluated. 2276 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. 2277 2278 .seealso: TSStep(), TSAdapt 2279 @*/ 2280 PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec U,PetscBool *done) 2281 { 2282 PetscErrorCode ierr; 2283 2284 PetscFunctionBegin; 2285 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2286 PetscValidType(ts,1); 2287 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 2288 if (!ts->ops->evaluatestep) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSEvaluateStep not implemented for type '%s'",((PetscObject)ts)->type_name); 2289 ierr = (*ts->ops->evaluatestep)(ts,order,U,done);CHKERRQ(ierr); 2290 PetscFunctionReturn(0); 2291 } 2292 2293 #undef __FUNCT__ 2294 #define __FUNCT__ "TSSolve" 2295 /*@ 2296 TSSolve - Steps the requested number of timesteps. 2297 2298 Collective on TS 2299 2300 Input Parameter: 2301 + ts - the TS context obtained from TSCreate() 2302 - u - the solution vector (can be null if TSSetSolution() was used, otherwise must contain the initial conditions) 2303 2304 Level: beginner 2305 2306 Notes: 2307 The final time returned by this function may be different from the time of the internally 2308 held state accessible by TSGetSolution() and TSGetTime() because the method may have 2309 stepped over the final time. 2310 2311 .keywords: TS, timestep, solve 2312 2313 .seealso: TSCreate(), TSSetSolution(), TSStep() 2314 @*/ 2315 PetscErrorCode TSSolve(TS ts,Vec u) 2316 { 2317 PetscBool flg; 2318 char filename[PETSC_MAX_PATH_LEN]; 2319 PetscViewer viewer; 2320 PetscErrorCode ierr; 2321 2322 PetscFunctionBegin; 2323 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2324 if (u) PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2325 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 */ 2326 if (!ts->vec_sol || u == ts->vec_sol) { 2327 Vec y; 2328 ierr = VecDuplicate(u,&y);CHKERRQ(ierr); 2329 ierr = TSSetSolution(ts,y);CHKERRQ(ierr); 2330 ierr = VecDestroy(&y);CHKERRQ(ierr); /* grant ownership */ 2331 } 2332 if (u) { 2333 ierr = VecCopy(u,ts->vec_sol);CHKERRQ(ierr); 2334 } 2335 } else { 2336 if (u) { 2337 ierr = TSSetSolution(ts,u);CHKERRQ(ierr); 2338 } 2339 } 2340 ierr = TSSetUp(ts);CHKERRQ(ierr); 2341 /* reset time step and iteration counters */ 2342 ts->steps = 0; 2343 ts->ksp_its = 0; 2344 ts->snes_its = 0; 2345 ts->num_snes_failures = 0; 2346 ts->reject = 0; 2347 ts->reason = TS_CONVERGED_ITERATING; 2348 2349 if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */ 2350 ierr = (*ts->ops->solve)(ts);CHKERRQ(ierr); 2351 ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr); 2352 ts->solvetime = ts->ptime; 2353 } else { 2354 /* steps the requested number of timesteps. */ 2355 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 2356 if (ts->steps >= ts->max_steps) 2357 ts->reason = TS_CONVERGED_ITS; 2358 else if (ts->ptime >= ts->max_time) 2359 ts->reason = TS_CONVERGED_TIME; 2360 while (!ts->reason) { 2361 ierr = TSStep(ts);CHKERRQ(ierr); 2362 ierr = TSPostStep(ts);CHKERRQ(ierr); 2363 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 2364 } 2365 if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) { 2366 ierr = TSInterpolate(ts,ts->max_time,u);CHKERRQ(ierr); 2367 ts->solvetime = ts->max_time; 2368 } else { 2369 ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr); 2370 ts->solvetime = ts->ptime; 2371 } 2372 } 2373 ierr = TSMonitor(ts,-1,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 2374 ierr = PetscOptionsGetString(((PetscObject)ts)->prefix,"-ts_view",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 2375 if (flg && !PetscPreLoadingOn) { 2376 ierr = PetscViewerASCIIOpen(((PetscObject)ts)->comm,filename,&viewer);CHKERRQ(ierr); 2377 ierr = TSView(ts,viewer);CHKERRQ(ierr); 2378 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2379 } 2380 flg = PETSC_FALSE; 2381 ierr = PetscOptionsGetBool(((PetscObject)ts)->prefix,"-ts_view_draw",&flg,PETSC_NULL);CHKERRQ(ierr); 2382 if (flg) { 2383 ierr = PetscViewerDrawOpen(((PetscObject)ts)->comm,PETSC_NULL,"TS Solver",0,0,600,600,&viewer);CHKERRQ(ierr); 2384 ierr = TSView(ts,viewer);CHKERRQ(ierr); 2385 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2386 } 2387 PetscFunctionReturn(0); 2388 } 2389 2390 #undef __FUNCT__ 2391 #define __FUNCT__ "TSMonitor" 2392 /*@ 2393 TSMonitor - Runs all user-provided monitor routines set using TSMonitorSet() 2394 2395 Collective on TS 2396 2397 Input Parameters: 2398 + ts - time stepping context obtained from TSCreate() 2399 . step - step number that has just completed 2400 . ptime - model time of the state 2401 - u - state at the current model time 2402 2403 Notes: 2404 TSMonitor() is typically used within the time stepping implementations. 2405 Users might call this function when using the TSStep() interface instead of TSSolve(). 2406 2407 Level: advanced 2408 2409 .keywords: TS, timestep 2410 @*/ 2411 PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u) 2412 { 2413 PetscErrorCode ierr; 2414 PetscInt i,n = ts->numbermonitors; 2415 2416 PetscFunctionBegin; 2417 for (i=0; i<n; i++) { 2418 ierr = (*ts->monitor[i])(ts,step,ptime,u,ts->monitorcontext[i]);CHKERRQ(ierr); 2419 } 2420 PetscFunctionReturn(0); 2421 } 2422 2423 /* ------------------------------------------------------------------------*/ 2424 struct _n_TSMonitorLGCtx { 2425 PetscDrawLG lg; 2426 PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */ 2427 PetscInt ksp_its,snes_its; 2428 }; 2429 2430 2431 #undef __FUNCT__ 2432 #define __FUNCT__ "TSMonitorLGCtxCreate" 2433 /*@C 2434 TSMonitorLGCtxCreate - Creates a line graph context for use with 2435 TS to monitor the solution process graphically in various ways 2436 2437 Collective on TS 2438 2439 Input Parameters: 2440 + host - the X display to open, or null for the local machine 2441 . label - the title to put in the title bar 2442 . x, y - the screen coordinates of the upper left coordinate of the window 2443 . m, n - the screen width and height in pixels 2444 - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time 2445 2446 Output Parameter: 2447 . ctx - the context 2448 2449 Options Database Key: 2450 + -ts_monitor_lg_timestep - automatically sets line graph monitor 2451 . -ts_monitor_lg_solution - 2452 . -ts_monitor_lg_error - 2453 . -ts_monitor_lg_ksp_iterations - 2454 . -ts_monitor_lg_snes_iterations - 2455 - -lg_indicate_data_points <true,false> - indicate the data points (at each time step) on the plot; default is true 2456 2457 Notes: 2458 Use TSMonitorLGCtxDestroy() to destroy. 2459 2460 Level: intermediate 2461 2462 .keywords: TS, monitor, line graph, residual, seealso 2463 2464 .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError() 2465 2466 @*/ 2467 PetscErrorCode TSMonitorLGCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorLGCtx *ctx) 2468 { 2469 PetscDraw win; 2470 PetscErrorCode ierr; 2471 PetscBool flg = PETSC_TRUE; 2472 2473 PetscFunctionBegin; 2474 ierr = PetscNew(struct _n_TSMonitorLGCtx,ctx);CHKERRQ(ierr); 2475 ierr = PetscDrawCreate(comm,host,label,x,y,m,n,&win);CHKERRQ(ierr); 2476 ierr = PetscDrawSetFromOptions(win);CHKERRQ(ierr); 2477 ierr = PetscDrawLGCreate(win,1,&(*ctx)->lg);CHKERRQ(ierr); 2478 ierr = PetscOptionsGetBool(PETSC_NULL,"-lg_indicate_data_points",&flg,PETSC_NULL);CHKERRQ(ierr); 2479 if (flg) { 2480 ierr = PetscDrawLGIndicateDataPoints((*ctx)->lg);CHKERRQ(ierr); 2481 } 2482 ierr = PetscLogObjectParent((*ctx)->lg,win);CHKERRQ(ierr); 2483 (*ctx)->howoften = howoften; 2484 PetscFunctionReturn(0); 2485 } 2486 2487 #undef __FUNCT__ 2488 #define __FUNCT__ "TSMonitorLGTimeStep" 2489 PetscErrorCode TSMonitorLGTimeStep(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 2490 { 2491 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 2492 PetscReal x = ptime,y; 2493 PetscErrorCode ierr; 2494 2495 PetscFunctionBegin; 2496 if (!n) { 2497 PetscDrawAxis axis; 2498 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 2499 ierr = PetscDrawAxisSetLabels(axis,"Timestep as function of time","Time","Time step");CHKERRQ(ierr); 2500 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 2501 } 2502 ierr = TSGetTimeStep(ts,&y);CHKERRQ(ierr); 2503 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 2504 if (((ctx->howoften > 0) && (!(n % ctx->howoften))) || ((ctx->howoften == -1) && (n == -1))){ 2505 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 2506 } 2507 PetscFunctionReturn(0); 2508 } 2509 2510 #undef __FUNCT__ 2511 #define __FUNCT__ "TSMonitorLGCtxDestroy" 2512 /*@C 2513 TSMonitorLGCtxDestroy - Destroys a line graph context that was created 2514 with TSMonitorLGCtxCreate(). 2515 2516 Collective on TSMonitorLGCtx 2517 2518 Input Parameter: 2519 . ctx - the monitor context 2520 2521 Level: intermediate 2522 2523 .keywords: TS, monitor, line graph, destroy 2524 2525 .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep(); 2526 @*/ 2527 PetscErrorCode TSMonitorLGCtxDestroy(TSMonitorLGCtx *ctx) 2528 { 2529 PetscDraw draw; 2530 PetscErrorCode ierr; 2531 2532 PetscFunctionBegin; 2533 ierr = PetscDrawLGGetDraw((*ctx)->lg,&draw);CHKERRQ(ierr); 2534 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 2535 ierr = PetscDrawLGDestroy(&(*ctx)->lg);CHKERRQ(ierr); 2536 ierr = PetscFree(*ctx);CHKERRQ(ierr); 2537 PetscFunctionReturn(0); 2538 } 2539 2540 #undef __FUNCT__ 2541 #define __FUNCT__ "TSGetTime" 2542 /*@ 2543 TSGetTime - Gets the time of the most recently completed step. 2544 2545 Not Collective 2546 2547 Input Parameter: 2548 . ts - the TS context obtained from TSCreate() 2549 2550 Output Parameter: 2551 . t - the current time 2552 2553 Level: beginner 2554 2555 Note: 2556 When called during time step evaluation (e.g. during residual evaluation or via hooks set using TSSetPreStep(), 2557 TSSetPreStage(), or TSSetPostStep()), the time is the time at the start of the step being evaluated. 2558 2559 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 2560 2561 .keywords: TS, get, time 2562 @*/ 2563 PetscErrorCode TSGetTime(TS ts,PetscReal* t) 2564 { 2565 PetscFunctionBegin; 2566 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2567 PetscValidRealPointer(t,2); 2568 *t = ts->ptime; 2569 PetscFunctionReturn(0); 2570 } 2571 2572 #undef __FUNCT__ 2573 #define __FUNCT__ "TSSetTime" 2574 /*@ 2575 TSSetTime - Allows one to reset the time. 2576 2577 Logically Collective on TS 2578 2579 Input Parameters: 2580 + ts - the TS context obtained from TSCreate() 2581 - time - the time 2582 2583 Level: intermediate 2584 2585 .seealso: TSGetTime(), TSSetDuration() 2586 2587 .keywords: TS, set, time 2588 @*/ 2589 PetscErrorCode TSSetTime(TS ts, PetscReal t) 2590 { 2591 PetscFunctionBegin; 2592 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2593 PetscValidLogicalCollectiveReal(ts,t,2); 2594 ts->ptime = t; 2595 PetscFunctionReturn(0); 2596 } 2597 2598 #undef __FUNCT__ 2599 #define __FUNCT__ "TSSetOptionsPrefix" 2600 /*@C 2601 TSSetOptionsPrefix - Sets the prefix used for searching for all 2602 TS options in the database. 2603 2604 Logically Collective on TS 2605 2606 Input Parameter: 2607 + ts - The TS context 2608 - prefix - The prefix to prepend to all option names 2609 2610 Notes: 2611 A hyphen (-) must NOT be given at the beginning of the prefix name. 2612 The first character of all runtime options is AUTOMATICALLY the 2613 hyphen. 2614 2615 Level: advanced 2616 2617 .keywords: TS, set, options, prefix, database 2618 2619 .seealso: TSSetFromOptions() 2620 2621 @*/ 2622 PetscErrorCode TSSetOptionsPrefix(TS ts,const char prefix[]) 2623 { 2624 PetscErrorCode ierr; 2625 SNES snes; 2626 2627 PetscFunctionBegin; 2628 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2629 ierr = PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2630 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2631 ierr = SNESSetOptionsPrefix(snes,prefix);CHKERRQ(ierr); 2632 PetscFunctionReturn(0); 2633 } 2634 2635 2636 #undef __FUNCT__ 2637 #define __FUNCT__ "TSAppendOptionsPrefix" 2638 /*@C 2639 TSAppendOptionsPrefix - Appends to the prefix used for searching for all 2640 TS options in the database. 2641 2642 Logically Collective on TS 2643 2644 Input Parameter: 2645 + ts - The TS context 2646 - prefix - The prefix to prepend to all option names 2647 2648 Notes: 2649 A hyphen (-) must NOT be given at the beginning of the prefix name. 2650 The first character of all runtime options is AUTOMATICALLY the 2651 hyphen. 2652 2653 Level: advanced 2654 2655 .keywords: TS, append, options, prefix, database 2656 2657 .seealso: TSGetOptionsPrefix() 2658 2659 @*/ 2660 PetscErrorCode TSAppendOptionsPrefix(TS ts,const char prefix[]) 2661 { 2662 PetscErrorCode ierr; 2663 SNES snes; 2664 2665 PetscFunctionBegin; 2666 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2667 ierr = PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2668 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2669 ierr = SNESAppendOptionsPrefix(snes,prefix);CHKERRQ(ierr); 2670 PetscFunctionReturn(0); 2671 } 2672 2673 #undef __FUNCT__ 2674 #define __FUNCT__ "TSGetOptionsPrefix" 2675 /*@C 2676 TSGetOptionsPrefix - Sets the prefix used for searching for all 2677 TS options in the database. 2678 2679 Not Collective 2680 2681 Input Parameter: 2682 . ts - The TS context 2683 2684 Output Parameter: 2685 . prefix - A pointer to the prefix string used 2686 2687 Notes: On the fortran side, the user should pass in a string 'prifix' of 2688 sufficient length to hold the prefix. 2689 2690 Level: intermediate 2691 2692 .keywords: TS, get, options, prefix, database 2693 2694 .seealso: TSAppendOptionsPrefix() 2695 @*/ 2696 PetscErrorCode TSGetOptionsPrefix(TS ts,const char *prefix[]) 2697 { 2698 PetscErrorCode ierr; 2699 2700 PetscFunctionBegin; 2701 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2702 PetscValidPointer(prefix,2); 2703 ierr = PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2704 PetscFunctionReturn(0); 2705 } 2706 2707 #undef __FUNCT__ 2708 #define __FUNCT__ "TSGetRHSJacobian" 2709 /*@C 2710 TSGetRHSJacobian - Returns the Jacobian J at the present timestep. 2711 2712 Not Collective, but parallel objects are returned if TS is parallel 2713 2714 Input Parameter: 2715 . ts - The TS context obtained from TSCreate() 2716 2717 Output Parameters: 2718 + J - The Jacobian J of F, where U_t = G(U,t) 2719 . M - The preconditioner matrix, usually the same as J 2720 . func - Function to compute the Jacobian of the RHS 2721 - ctx - User-defined context for Jacobian evaluation routine 2722 2723 Notes: You can pass in PETSC_NULL for any return argument you do not need. 2724 2725 Level: intermediate 2726 2727 .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 2728 2729 .keywords: TS, timestep, get, matrix, Jacobian 2730 @*/ 2731 PetscErrorCode TSGetRHSJacobian(TS ts,Mat *J,Mat *M,TSRHSJacobian *func,void **ctx) 2732 { 2733 PetscErrorCode ierr; 2734 SNES snes; 2735 DM dm; 2736 2737 PetscFunctionBegin; 2738 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2739 ierr = SNESGetJacobian(snes,J,M,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2740 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 2741 ierr = DMTSGetRHSJacobian(dm,func,ctx);CHKERRQ(ierr); 2742 PetscFunctionReturn(0); 2743 } 2744 2745 #undef __FUNCT__ 2746 #define __FUNCT__ "TSGetIJacobian" 2747 /*@C 2748 TSGetIJacobian - Returns the implicit Jacobian at the present timestep. 2749 2750 Not Collective, but parallel objects are returned if TS is parallel 2751 2752 Input Parameter: 2753 . ts - The TS context obtained from TSCreate() 2754 2755 Output Parameters: 2756 + A - The Jacobian of F(t,U,U_t) 2757 . B - The preconditioner matrix, often the same as A 2758 . f - The function to compute the matrices 2759 - ctx - User-defined context for Jacobian evaluation routine 2760 2761 Notes: You can pass in PETSC_NULL for any return argument you do not need. 2762 2763 Level: advanced 2764 2765 .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 2766 2767 .keywords: TS, timestep, get, matrix, Jacobian 2768 @*/ 2769 PetscErrorCode TSGetIJacobian(TS ts,Mat *A,Mat *B,TSIJacobian *f,void **ctx) 2770 { 2771 PetscErrorCode ierr; 2772 SNES snes; 2773 DM dm; 2774 2775 PetscFunctionBegin; 2776 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2777 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 2778 ierr = SNESGetJacobian(snes,A,B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2779 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 2780 ierr = DMTSGetIJacobian(dm,f,ctx);CHKERRQ(ierr); 2781 PetscFunctionReturn(0); 2782 } 2783 2784 struct _n_TSMonitorDrawCtx { 2785 PetscViewer viewer; 2786 Vec initialsolution; 2787 PetscBool showinitial; 2788 PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */ 2789 }; 2790 2791 #undef __FUNCT__ 2792 #define __FUNCT__ "TSMonitorDrawSolution" 2793 /*@C 2794 TSMonitorDrawSolution - Monitors progress of the TS solvers by calling 2795 VecView() for the solution at each timestep 2796 2797 Collective on TS 2798 2799 Input Parameters: 2800 + ts - the TS context 2801 . step - current time-step 2802 . ptime - current time 2803 - dummy - either a viewer or PETSC_NULL 2804 2805 Options Database: 2806 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 2807 2808 Notes: the initial solution and current solution are not displayed with a common axis scaling so generally the option -ts_monitor_draw_solution_initial 2809 will look bad 2810 2811 Level: intermediate 2812 2813 .keywords: TS, vector, monitor, view 2814 2815 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 2816 @*/ 2817 PetscErrorCode TSMonitorDrawSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 2818 { 2819 PetscErrorCode ierr; 2820 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy; 2821 2822 PetscFunctionBegin; 2823 if (!step && ictx->showinitial) { 2824 if (!ictx->initialsolution) { 2825 ierr = VecDuplicate(u,&ictx->initialsolution);CHKERRQ(ierr); 2826 } 2827 ierr = VecCopy(u,ictx->initialsolution);CHKERRQ(ierr); 2828 } 2829 if (!(((ictx->howoften > 0) && (!(step % ictx->howoften)) && (step > -1)) || ((ictx->howoften == -1) && (step == -1)))) PetscFunctionReturn(0); 2830 2831 if (ictx->showinitial) { 2832 PetscReal pause; 2833 ierr = PetscViewerDrawGetPause(ictx->viewer,&pause);CHKERRQ(ierr); 2834 ierr = PetscViewerDrawSetPause(ictx->viewer,0.0);CHKERRQ(ierr); 2835 ierr = VecView(ictx->initialsolution,ictx->viewer);CHKERRQ(ierr); 2836 ierr = PetscViewerDrawSetPause(ictx->viewer,pause);CHKERRQ(ierr); 2837 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);CHKERRQ(ierr); 2838 } 2839 ierr = VecView(u,ictx->viewer);CHKERRQ(ierr); 2840 if (ictx->showinitial) { 2841 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);CHKERRQ(ierr); 2842 } 2843 PetscFunctionReturn(0); 2844 } 2845 2846 2847 #undef __FUNCT__ 2848 #define __FUNCT__ "TSMonitorDrawCtxDestroy" 2849 /*@C 2850 TSMonitorDrawCtxDestroy - Destroys the monitor context for TSMonitorDrawSolution() 2851 2852 Collective on TS 2853 2854 Input Parameters: 2855 . ctx - the monitor context 2856 2857 Level: intermediate 2858 2859 .keywords: TS, vector, monitor, view 2860 2861 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawSolution(), TSMonitorDrawError() 2862 @*/ 2863 PetscErrorCode TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx) 2864 { 2865 PetscErrorCode ierr; 2866 2867 PetscFunctionBegin; 2868 ierr = PetscViewerDestroy(&(*ictx)->viewer);CHKERRQ(ierr); 2869 ierr = VecDestroy(&(*ictx)->initialsolution);CHKERRQ(ierr); 2870 ierr = PetscFree(*ictx);CHKERRQ(ierr); 2871 PetscFunctionReturn(0); 2872 } 2873 2874 #undef __FUNCT__ 2875 #define __FUNCT__ "TSMonitorDrawCtxCreate" 2876 /*@C 2877 TSMonitorDrawCtxCreate - Creates the monitor context for TSMonitorDrawCtx 2878 2879 Collective on TS 2880 2881 Input Parameter: 2882 . ts - time-step context 2883 2884 Output Patameter: 2885 . ctx - the monitor context 2886 2887 Options Database: 2888 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 2889 2890 Level: intermediate 2891 2892 .keywords: TS, vector, monitor, view 2893 2894 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawCtx() 2895 @*/ 2896 PetscErrorCode TSMonitorDrawCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorDrawCtx *ctx) 2897 { 2898 PetscErrorCode ierr; 2899 2900 PetscFunctionBegin; 2901 ierr = PetscNew(struct _n_TSMonitorDrawCtx,ctx);CHKERRQ(ierr); 2902 ierr = PetscViewerDrawOpen(comm,host,label,x,y,m,n,&(*ctx)->viewer);CHKERRQ(ierr); 2903 (*ctx)->showinitial = PETSC_FALSE; 2904 (*ctx)->howoften = howoften; 2905 ierr = PetscOptionsGetBool(PETSC_NULL,"-ts_monitor_draw_solution_initial",&(*ctx)->showinitial,PETSC_NULL);CHKERRQ(ierr); 2906 PetscFunctionReturn(0); 2907 } 2908 2909 #undef __FUNCT__ 2910 #define __FUNCT__ "TSMonitorDrawError" 2911 /*@C 2912 TSMonitorDrawError - Monitors progress of the TS solvers by calling 2913 VecView() for the error at each timestep 2914 2915 Collective on TS 2916 2917 Input Parameters: 2918 + ts - the TS context 2919 . step - current time-step 2920 . ptime - current time 2921 - dummy - either a viewer or PETSC_NULL 2922 2923 Level: intermediate 2924 2925 .keywords: TS, vector, monitor, view 2926 2927 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 2928 @*/ 2929 PetscErrorCode TSMonitorDrawError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 2930 { 2931 PetscErrorCode ierr; 2932 TSMonitorDrawCtx ctx = (TSMonitorDrawCtx)dummy; 2933 PetscViewer viewer = ctx->viewer; 2934 Vec work; 2935 2936 PetscFunctionBegin; 2937 if (!(((ctx->howoften > 0) && (!(step % ctx->howoften)) && (step > -1)) || ((ctx->howoften == -1) && (step == -1)))) PetscFunctionReturn(0); 2938 ierr = VecDuplicate(u,&work);CHKERRQ(ierr); 2939 ierr = TSComputeSolutionFunction(ts,ptime,work);CHKERRQ(ierr); 2940 ierr = VecAXPY(work,-1.0,u);CHKERRQ(ierr); 2941 ierr = VecView(work,viewer);CHKERRQ(ierr); 2942 ierr = VecDestroy(&work);CHKERRQ(ierr); 2943 PetscFunctionReturn(0); 2944 } 2945 2946 #include <petsc-private/dmimpl.h> 2947 #undef __FUNCT__ 2948 #define __FUNCT__ "TSSetDM" 2949 /*@ 2950 TSSetDM - Sets the DM that may be used by some preconditioners 2951 2952 Logically Collective on TS and DM 2953 2954 Input Parameters: 2955 + ts - the preconditioner context 2956 - dm - the dm 2957 2958 Level: intermediate 2959 2960 2961 .seealso: TSGetDM(), SNESSetDM(), SNESGetDM() 2962 @*/ 2963 PetscErrorCode TSSetDM(TS ts,DM dm) 2964 { 2965 PetscErrorCode ierr; 2966 SNES snes; 2967 DMTS tsdm; 2968 2969 PetscFunctionBegin; 2970 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2971 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 2972 if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */ 2973 if (ts->dm->dmts && !dm->dmts) { 2974 ierr = DMCopyDMTS(ts->dm,dm);CHKERRQ(ierr); 2975 ierr = DMGetDMTS(ts->dm,&tsdm);CHKERRQ(ierr); 2976 if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */ 2977 tsdm->originaldm = dm; 2978 } 2979 } 2980 ierr = DMDestroy(&ts->dm);CHKERRQ(ierr); 2981 } 2982 ts->dm = dm; 2983 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2984 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 2985 PetscFunctionReturn(0); 2986 } 2987 2988 #undef __FUNCT__ 2989 #define __FUNCT__ "TSGetDM" 2990 /*@ 2991 TSGetDM - Gets the DM that may be used by some preconditioners 2992 2993 Not Collective 2994 2995 Input Parameter: 2996 . ts - the preconditioner context 2997 2998 Output Parameter: 2999 . dm - the dm 3000 3001 Level: intermediate 3002 3003 3004 .seealso: TSSetDM(), SNESSetDM(), SNESGetDM() 3005 @*/ 3006 PetscErrorCode TSGetDM(TS ts,DM *dm) 3007 { 3008 PetscErrorCode ierr; 3009 3010 PetscFunctionBegin; 3011 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3012 if (!ts->dm) { 3013 ierr = DMShellCreate(((PetscObject)ts)->comm,&ts->dm);CHKERRQ(ierr); 3014 if (ts->snes) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 3015 } 3016 *dm = ts->dm; 3017 PetscFunctionReturn(0); 3018 } 3019 3020 #undef __FUNCT__ 3021 #define __FUNCT__ "SNESTSFormFunction" 3022 /*@ 3023 SNESTSFormFunction - Function to evaluate nonlinear residual 3024 3025 Logically Collective on SNES 3026 3027 Input Parameter: 3028 + snes - nonlinear solver 3029 . U - the current state at which to evaluate the residual 3030 - ctx - user context, must be a TS 3031 3032 Output Parameter: 3033 . F - the nonlinear residual 3034 3035 Notes: 3036 This function is not normally called by users and is automatically registered with the SNES used by TS. 3037 It is most frequently passed to MatFDColoringSetFunction(). 3038 3039 Level: advanced 3040 3041 .seealso: SNESSetFunction(), MatFDColoringSetFunction() 3042 @*/ 3043 PetscErrorCode SNESTSFormFunction(SNES snes,Vec U,Vec F,void *ctx) 3044 { 3045 TS ts = (TS)ctx; 3046 PetscErrorCode ierr; 3047 3048 PetscFunctionBegin; 3049 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3050 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 3051 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 3052 PetscValidHeaderSpecific(ts,TS_CLASSID,4); 3053 ierr = (ts->ops->snesfunction)(snes,U,F,ts);CHKERRQ(ierr); 3054 PetscFunctionReturn(0); 3055 } 3056 3057 #undef __FUNCT__ 3058 #define __FUNCT__ "SNESTSFormJacobian" 3059 /*@ 3060 SNESTSFormJacobian - Function to evaluate the Jacobian 3061 3062 Collective on SNES 3063 3064 Input Parameter: 3065 + snes - nonlinear solver 3066 . U - the current state at which to evaluate the residual 3067 - ctx - user context, must be a TS 3068 3069 Output Parameter: 3070 + A - the Jacobian 3071 . B - the preconditioning matrix (may be the same as A) 3072 - flag - indicates any structure change in the matrix 3073 3074 Notes: 3075 This function is not normally called by users and is automatically registered with the SNES used by TS. 3076 3077 Level: developer 3078 3079 .seealso: SNESSetJacobian() 3080 @*/ 3081 PetscErrorCode SNESTSFormJacobian(SNES snes,Vec U,Mat *A,Mat *B,MatStructure *flag,void *ctx) 3082 { 3083 TS ts = (TS)ctx; 3084 PetscErrorCode ierr; 3085 3086 PetscFunctionBegin; 3087 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3088 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 3089 PetscValidPointer(A,3); 3090 PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 3091 PetscValidPointer(B,4); 3092 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); 3093 PetscValidPointer(flag,5); 3094 PetscValidHeaderSpecific(ts,TS_CLASSID,6); 3095 ierr = (ts->ops->snesjacobian)(snes,U,A,B,flag,ts);CHKERRQ(ierr); 3096 PetscFunctionReturn(0); 3097 } 3098 3099 #undef __FUNCT__ 3100 #define __FUNCT__ "TSComputeRHSFunctionLinear" 3101 /*@C 3102 TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems only 3103 3104 Collective on TS 3105 3106 Input Arguments: 3107 + ts - time stepping context 3108 . t - time at which to evaluate 3109 . U - state at which to evaluate 3110 - ctx - context 3111 3112 Output Arguments: 3113 . F - right hand side 3114 3115 Level: intermediate 3116 3117 Notes: 3118 This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems. 3119 The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian(). 3120 3121 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant() 3122 @*/ 3123 PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec U,Vec F,void *ctx) 3124 { 3125 PetscErrorCode ierr; 3126 Mat Arhs,Brhs; 3127 MatStructure flg2; 3128 3129 PetscFunctionBegin; 3130 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 3131 ierr = TSComputeRHSJacobian(ts,t,U,&Arhs,&Brhs,&flg2);CHKERRQ(ierr); 3132 ierr = MatMult(Arhs,U,F);CHKERRQ(ierr); 3133 PetscFunctionReturn(0); 3134 } 3135 3136 #undef __FUNCT__ 3137 #define __FUNCT__ "TSComputeRHSJacobianConstant" 3138 /*@C 3139 TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent. 3140 3141 Collective on TS 3142 3143 Input Arguments: 3144 + ts - time stepping context 3145 . t - time at which to evaluate 3146 . U - state at which to evaluate 3147 - ctx - context 3148 3149 Output Arguments: 3150 + A - pointer to operator 3151 . B - pointer to preconditioning matrix 3152 - flg - matrix structure flag 3153 3154 Level: intermediate 3155 3156 Notes: 3157 This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems. 3158 3159 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear() 3160 @*/ 3161 PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec U,Mat *A,Mat *B,MatStructure *flg,void *ctx) 3162 { 3163 PetscFunctionBegin; 3164 *flg = SAME_PRECONDITIONER; 3165 PetscFunctionReturn(0); 3166 } 3167 3168 #undef __FUNCT__ 3169 #define __FUNCT__ "TSComputeIFunctionLinear" 3170 /*@C 3171 TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only 3172 3173 Collective on TS 3174 3175 Input Arguments: 3176 + ts - time stepping context 3177 . t - time at which to evaluate 3178 . U - state at which to evaluate 3179 . Udot - time derivative of state vector 3180 - ctx - context 3181 3182 Output Arguments: 3183 . F - left hand side 3184 3185 Level: intermediate 3186 3187 Notes: 3188 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 3189 user is required to write their own TSComputeIFunction. 3190 This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems. 3191 The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian(). 3192 3193 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant() 3194 @*/ 3195 PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx) 3196 { 3197 PetscErrorCode ierr; 3198 Mat A,B; 3199 MatStructure flg2; 3200 3201 PetscFunctionBegin; 3202 ierr = TSGetIJacobian(ts,&A,&B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 3203 ierr = TSComputeIJacobian(ts,t,U,Udot,1.0,&A,&B,&flg2,PETSC_TRUE);CHKERRQ(ierr); 3204 ierr = MatMult(A,Udot,F);CHKERRQ(ierr); 3205 PetscFunctionReturn(0); 3206 } 3207 3208 #undef __FUNCT__ 3209 #define __FUNCT__ "TSComputeIJacobianConstant" 3210 /*@C 3211 TSComputeIJacobianConstant - Reuses a Jacobian that is time-independent. 3212 3213 Collective on TS 3214 3215 Input Arguments: 3216 + ts - time stepping context 3217 . t - time at which to evaluate 3218 . U - state at which to evaluate 3219 . Udot - time derivative of state vector 3220 . shift - shift to apply 3221 - ctx - context 3222 3223 Output Arguments: 3224 + A - pointer to operator 3225 . B - pointer to preconditioning matrix 3226 - flg - matrix structure flag 3227 3228 Level: intermediate 3229 3230 Notes: 3231 This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems. 3232 3233 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear() 3234 @*/ 3235 PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,void *ctx) 3236 { 3237 PetscFunctionBegin; 3238 *flg = SAME_PRECONDITIONER; 3239 PetscFunctionReturn(0); 3240 } 3241 3242 #undef __FUNCT__ 3243 #define __FUNCT__ "TSGetConvergedReason" 3244 /*@ 3245 TSGetConvergedReason - Gets the reason the TS iteration was stopped. 3246 3247 Not Collective 3248 3249 Input Parameter: 3250 . ts - the TS context 3251 3252 Output Parameter: 3253 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 3254 manual pages for the individual convergence tests for complete lists 3255 3256 Level: beginner 3257 3258 Notes: 3259 Can only be called after the call to TSSolve() is complete. 3260 3261 .keywords: TS, nonlinear, set, convergence, test 3262 3263 .seealso: TSSetConvergenceTest(), TSConvergedReason 3264 @*/ 3265 PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason) 3266 { 3267 PetscFunctionBegin; 3268 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3269 PetscValidPointer(reason,2); 3270 *reason = ts->reason; 3271 PetscFunctionReturn(0); 3272 } 3273 3274 #undef __FUNCT__ 3275 #define __FUNCT__ "TSGetSolveTime" 3276 /*@ 3277 TSGetSolveTime - Gets the time after a call to TSSolve() 3278 3279 Not Collective 3280 3281 Input Parameter: 3282 . ts - the TS context 3283 3284 Output Parameter: 3285 . ftime - the final time. This time should correspond to the final time set with TSSetDuration() 3286 3287 Level: beginner 3288 3289 Notes: 3290 Can only be called after the call to TSSolve() is complete. 3291 3292 .keywords: TS, nonlinear, set, convergence, test 3293 3294 .seealso: TSSetConvergenceTest(), TSConvergedReason 3295 @*/ 3296 PetscErrorCode TSGetSolveTime(TS ts,PetscReal *ftime) 3297 { 3298 PetscFunctionBegin; 3299 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3300 PetscValidPointer(ftime,2); 3301 *ftime = ts->solvetime; 3302 PetscFunctionReturn(0); 3303 } 3304 3305 #undef __FUNCT__ 3306 #define __FUNCT__ "TSGetSNESIterations" 3307 /*@ 3308 TSGetSNESIterations - Gets the total number of nonlinear iterations 3309 used by the time integrator. 3310 3311 Not Collective 3312 3313 Input Parameter: 3314 . ts - TS context 3315 3316 Output Parameter: 3317 . nits - number of nonlinear iterations 3318 3319 Notes: 3320 This counter is reset to zero for each successive call to TSSolve(). 3321 3322 Level: intermediate 3323 3324 .keywords: TS, get, number, nonlinear, iterations 3325 3326 .seealso: TSGetKSPIterations() 3327 @*/ 3328 PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits) 3329 { 3330 PetscFunctionBegin; 3331 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3332 PetscValidIntPointer(nits,2); 3333 *nits = ts->snes_its; 3334 PetscFunctionReturn(0); 3335 } 3336 3337 #undef __FUNCT__ 3338 #define __FUNCT__ "TSGetKSPIterations" 3339 /*@ 3340 TSGetKSPIterations - Gets the total number of linear iterations 3341 used by the time integrator. 3342 3343 Not Collective 3344 3345 Input Parameter: 3346 . ts - TS context 3347 3348 Output Parameter: 3349 . lits - number of linear iterations 3350 3351 Notes: 3352 This counter is reset to zero for each successive call to TSSolve(). 3353 3354 Level: intermediate 3355 3356 .keywords: TS, get, number, linear, iterations 3357 3358 .seealso: TSGetSNESIterations(), SNESGetKSPIterations() 3359 @*/ 3360 PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits) 3361 { 3362 PetscFunctionBegin; 3363 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3364 PetscValidIntPointer(lits,2); 3365 *lits = ts->ksp_its; 3366 PetscFunctionReturn(0); 3367 } 3368 3369 #undef __FUNCT__ 3370 #define __FUNCT__ "TSGetStepRejections" 3371 /*@ 3372 TSGetStepRejections - Gets the total number of rejected steps. 3373 3374 Not Collective 3375 3376 Input Parameter: 3377 . ts - TS context 3378 3379 Output Parameter: 3380 . rejects - number of steps rejected 3381 3382 Notes: 3383 This counter is reset to zero for each successive call to TSSolve(). 3384 3385 Level: intermediate 3386 3387 .keywords: TS, get, number 3388 3389 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails() 3390 @*/ 3391 PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects) 3392 { 3393 PetscFunctionBegin; 3394 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3395 PetscValidIntPointer(rejects,2); 3396 *rejects = ts->reject; 3397 PetscFunctionReturn(0); 3398 } 3399 3400 #undef __FUNCT__ 3401 #define __FUNCT__ "TSGetSNESFailures" 3402 /*@ 3403 TSGetSNESFailures - Gets the total number of failed SNES solves 3404 3405 Not Collective 3406 3407 Input Parameter: 3408 . ts - TS context 3409 3410 Output Parameter: 3411 . fails - number of failed nonlinear solves 3412 3413 Notes: 3414 This counter is reset to zero for each successive call to TSSolve(). 3415 3416 Level: intermediate 3417 3418 .keywords: TS, get, number 3419 3420 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures() 3421 @*/ 3422 PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails) 3423 { 3424 PetscFunctionBegin; 3425 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3426 PetscValidIntPointer(fails,2); 3427 *fails = ts->num_snes_failures; 3428 PetscFunctionReturn(0); 3429 } 3430 3431 #undef __FUNCT__ 3432 #define __FUNCT__ "TSSetMaxStepRejections" 3433 /*@ 3434 TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails 3435 3436 Not Collective 3437 3438 Input Parameter: 3439 + ts - TS context 3440 - rejects - maximum number of rejected steps, pass -1 for unlimited 3441 3442 Notes: 3443 The counter is reset to zero for each step 3444 3445 Options Database Key: 3446 . -ts_max_reject - Maximum number of step rejections before a step fails 3447 3448 Level: intermediate 3449 3450 .keywords: TS, set, maximum, number 3451 3452 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 3453 @*/ 3454 PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects) 3455 { 3456 PetscFunctionBegin; 3457 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3458 ts->max_reject = rejects; 3459 PetscFunctionReturn(0); 3460 } 3461 3462 #undef __FUNCT__ 3463 #define __FUNCT__ "TSSetMaxSNESFailures" 3464 /*@ 3465 TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves 3466 3467 Not Collective 3468 3469 Input Parameter: 3470 + ts - TS context 3471 - fails - maximum number of failed nonlinear solves, pass -1 for unlimited 3472 3473 Notes: 3474 The counter is reset to zero for each successive call to TSSolve(). 3475 3476 Options Database Key: 3477 . -ts_max_snes_failures - Maximum number of nonlinear solve failures 3478 3479 Level: intermediate 3480 3481 .keywords: TS, set, maximum, number 3482 3483 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason() 3484 @*/ 3485 PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails) 3486 { 3487 PetscFunctionBegin; 3488 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3489 ts->max_snes_failures = fails; 3490 PetscFunctionReturn(0); 3491 } 3492 3493 #undef __FUNCT__ 3494 #define __FUNCT__ "TSSetErrorIfStepFails()" 3495 /*@ 3496 TSSetErrorIfStepFails - Error if no step succeeds 3497 3498 Not Collective 3499 3500 Input Parameter: 3501 + ts - TS context 3502 - err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure 3503 3504 Options Database Key: 3505 . -ts_error_if_step_fails - Error if no step succeeds 3506 3507 Level: intermediate 3508 3509 .keywords: TS, set, error 3510 3511 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 3512 @*/ 3513 PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err) 3514 { 3515 PetscFunctionBegin; 3516 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3517 ts->errorifstepfailed = err; 3518 PetscFunctionReturn(0); 3519 } 3520 3521 #undef __FUNCT__ 3522 #define __FUNCT__ "TSMonitorSolutionBinary" 3523 /*@C 3524 TSMonitorSolutionBinary - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file 3525 3526 Collective on TS 3527 3528 Input Parameters: 3529 + ts - the TS context 3530 . step - current time-step 3531 . ptime - current time 3532 . u - current state 3533 - viewer - binary viewer 3534 3535 Level: intermediate 3536 3537 .keywords: TS, vector, monitor, view 3538 3539 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3540 @*/ 3541 PetscErrorCode TSMonitorSolutionBinary(TS ts,PetscInt step,PetscReal ptime,Vec u,void *viewer) 3542 { 3543 PetscErrorCode ierr; 3544 PetscViewer v = (PetscViewer)viewer; 3545 3546 PetscFunctionBegin; 3547 ierr = VecView(u,v);CHKERRQ(ierr); 3548 PetscFunctionReturn(0); 3549 } 3550 3551 #undef __FUNCT__ 3552 #define __FUNCT__ "TSMonitorSolutionVTK" 3553 /*@C 3554 TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep. 3555 3556 Collective on TS 3557 3558 Input Parameters: 3559 + ts - the TS context 3560 . step - current time-step 3561 . ptime - current time 3562 . u - current state 3563 - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 3564 3565 Level: intermediate 3566 3567 Notes: 3568 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. 3569 These are named according to the file name template. 3570 3571 This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy(). 3572 3573 .keywords: TS, vector, monitor, view 3574 3575 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3576 @*/ 3577 PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec u,void *filenametemplate) 3578 { 3579 PetscErrorCode ierr; 3580 char filename[PETSC_MAX_PATH_LEN]; 3581 PetscViewer viewer; 3582 3583 PetscFunctionBegin; 3584 ierr = PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);CHKERRQ(ierr); 3585 ierr = PetscViewerVTKOpen(((PetscObject)ts)->comm,filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 3586 ierr = VecView(u,viewer);CHKERRQ(ierr); 3587 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3588 PetscFunctionReturn(0); 3589 } 3590 3591 #undef __FUNCT__ 3592 #define __FUNCT__ "TSMonitorSolutionVTKDestroy" 3593 /*@C 3594 TSMonitorSolutionVTKDestroy - Destroy context for monitoring 3595 3596 Collective on TS 3597 3598 Input Parameters: 3599 . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 3600 3601 Level: intermediate 3602 3603 Note: 3604 This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK(). 3605 3606 .keywords: TS, vector, monitor, view 3607 3608 .seealso: TSMonitorSet(), TSMonitorSolutionVTK() 3609 @*/ 3610 PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate) 3611 { 3612 PetscErrorCode ierr; 3613 3614 PetscFunctionBegin; 3615 ierr = PetscFree(*(char**)filenametemplate);CHKERRQ(ierr); 3616 PetscFunctionReturn(0); 3617 } 3618 3619 #undef __FUNCT__ 3620 #define __FUNCT__ "TSGetAdapt" 3621 /*@ 3622 TSGetAdapt - Get the adaptive controller context for the current method 3623 3624 Collective on TS if controller has not been created yet 3625 3626 Input Arguments: 3627 . ts - time stepping context 3628 3629 Output Arguments: 3630 . adapt - adaptive controller 3631 3632 Level: intermediate 3633 3634 .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose() 3635 @*/ 3636 PetscErrorCode TSGetAdapt(TS ts,TSAdapt *adapt) 3637 { 3638 PetscErrorCode ierr; 3639 3640 PetscFunctionBegin; 3641 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3642 PetscValidPointer(adapt,2); 3643 if (!ts->adapt) { 3644 ierr = TSAdaptCreate(((PetscObject)ts)->comm,&ts->adapt);CHKERRQ(ierr); 3645 ierr = PetscLogObjectParent(ts,ts->adapt);CHKERRQ(ierr); 3646 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);CHKERRQ(ierr); 3647 } 3648 *adapt = ts->adapt; 3649 PetscFunctionReturn(0); 3650 } 3651 3652 #undef __FUNCT__ 3653 #define __FUNCT__ "TSSetTolerances" 3654 /*@ 3655 TSSetTolerances - Set tolerances for local truncation error when using adaptive controller 3656 3657 Logically Collective 3658 3659 Input Arguments: 3660 + ts - time integration context 3661 . atol - scalar absolute tolerances, PETSC_DECIDE to leave current value 3662 . vatol - vector of absolute tolerances or PETSC_NULL, used in preference to atol if present 3663 . rtol - scalar relative tolerances, PETSC_DECIDE to leave current value 3664 - vrtol - vector of relative tolerances or PETSC_NULL, used in preference to atol if present 3665 3666 Level: beginner 3667 3668 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSGetTolerances() 3669 @*/ 3670 PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol) 3671 { 3672 PetscErrorCode ierr; 3673 3674 PetscFunctionBegin; 3675 if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol; 3676 if (vatol) { 3677 ierr = PetscObjectReference((PetscObject)vatol);CHKERRQ(ierr); 3678 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 3679 ts->vatol = vatol; 3680 } 3681 if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol; 3682 if (vrtol) { 3683 ierr = PetscObjectReference((PetscObject)vrtol);CHKERRQ(ierr); 3684 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 3685 ts->vrtol = vrtol; 3686 } 3687 PetscFunctionReturn(0); 3688 } 3689 3690 #undef __FUNCT__ 3691 #define __FUNCT__ "TSGetTolerances" 3692 /*@ 3693 TSGetTolerances - Get tolerances for local truncation error when using adaptive controller 3694 3695 Logically Collective 3696 3697 Input Arguments: 3698 . ts - time integration context 3699 3700 Output Arguments: 3701 + atol - scalar absolute tolerances, PETSC_NULL to ignore 3702 . vatol - vector of absolute tolerances, PETSC_NULL to ignore 3703 . rtol - scalar relative tolerances, PETSC_NULL to ignore 3704 - vrtol - vector of relative tolerances, PETSC_NULL to ignore 3705 3706 Level: beginner 3707 3708 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSSetTolerances() 3709 @*/ 3710 PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol) 3711 { 3712 PetscFunctionBegin; 3713 if (atol) *atol = ts->atol; 3714 if (vatol) *vatol = ts->vatol; 3715 if (rtol) *rtol = ts->rtol; 3716 if (vrtol) *vrtol = ts->vrtol; 3717 PetscFunctionReturn(0); 3718 } 3719 3720 #undef __FUNCT__ 3721 #define __FUNCT__ "TSErrorNormWRMS" 3722 /*@ 3723 TSErrorNormWRMS - compute a weighted norm of the difference between a vector and the current state 3724 3725 Collective on TS 3726 3727 Input Arguments: 3728 + ts - time stepping context 3729 - Y - state vector to be compared to ts->vec_sol 3730 3731 Output Arguments: 3732 . norm - weighted norm, a value of 1.0 is considered small 3733 3734 Level: developer 3735 3736 .seealso: TSSetTolerances() 3737 @*/ 3738 PetscErrorCode TSErrorNormWRMS(TS ts,Vec Y,PetscReal *norm) 3739 { 3740 PetscErrorCode ierr; 3741 PetscInt i,n,N; 3742 const PetscScalar *u,*y; 3743 Vec U; 3744 PetscReal sum,gsum; 3745 3746 PetscFunctionBegin; 3747 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3748 PetscValidHeaderSpecific(Y,VEC_CLASSID,2); 3749 PetscValidPointer(norm,3); 3750 U = ts->vec_sol; 3751 PetscCheckSameTypeAndComm(U,1,Y,2); 3752 if (U == Y) SETERRQ(((PetscObject)U)->comm,PETSC_ERR_ARG_IDN,"Y cannot be the TS solution vector"); 3753 3754 ierr = VecGetSize(U,&N);CHKERRQ(ierr); 3755 ierr = VecGetLocalSize(U,&n);CHKERRQ(ierr); 3756 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 3757 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 3758 sum = 0.; 3759 if (ts->vatol && ts->vrtol) { 3760 const PetscScalar *atol,*rtol; 3761 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3762 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3763 for (i=0; i<n; i++) { 3764 PetscReal tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 3765 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 3766 } 3767 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3768 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3769 } else if (ts->vatol) { /* vector atol, scalar rtol */ 3770 const PetscScalar *atol; 3771 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3772 for (i=0; i<n; i++) { 3773 PetscReal tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 3774 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 3775 } 3776 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3777 } else if (ts->vrtol) { /* scalar atol, vector rtol */ 3778 const PetscScalar *rtol; 3779 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3780 for (i=0; i<n; i++) { 3781 PetscReal tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 3782 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 3783 } 3784 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3785 } else { /* scalar atol, scalar rtol */ 3786 for (i=0; i<n; i++) { 3787 PetscReal tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 3788 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 3789 } 3790 } 3791 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 3792 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 3793 3794 ierr = MPI_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,((PetscObject)ts)->comm);CHKERRQ(ierr); 3795 *norm = PetscSqrtReal(gsum / N); 3796 if (PetscIsInfOrNanScalar(*norm)) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 3797 PetscFunctionReturn(0); 3798 } 3799 3800 #undef __FUNCT__ 3801 #define __FUNCT__ "TSSetCFLTimeLocal" 3802 /*@ 3803 TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler 3804 3805 Logically Collective on TS 3806 3807 Input Arguments: 3808 + ts - time stepping context 3809 - cfltime - maximum stable time step if using forward Euler (value can be different on each process) 3810 3811 Note: 3812 After calling this function, the global CFL time can be obtained by calling TSGetCFLTime() 3813 3814 Level: intermediate 3815 3816 .seealso: TSGetCFLTime(), TSADAPTCFL 3817 @*/ 3818 PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime) 3819 { 3820 PetscFunctionBegin; 3821 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3822 ts->cfltime_local = cfltime; 3823 ts->cfltime = -1.; 3824 PetscFunctionReturn(0); 3825 } 3826 3827 #undef __FUNCT__ 3828 #define __FUNCT__ "TSGetCFLTime" 3829 /*@ 3830 TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler 3831 3832 Collective on TS 3833 3834 Input Arguments: 3835 . ts - time stepping context 3836 3837 Output Arguments: 3838 . cfltime - maximum stable time step for forward Euler 3839 3840 Level: advanced 3841 3842 .seealso: TSSetCFLTimeLocal() 3843 @*/ 3844 PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime) 3845 { 3846 PetscErrorCode ierr; 3847 3848 PetscFunctionBegin; 3849 if (ts->cfltime < 0) { 3850 ierr = MPI_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,((PetscObject)ts)->comm);CHKERRQ(ierr); 3851 } 3852 *cfltime = ts->cfltime; 3853 PetscFunctionReturn(0); 3854 } 3855 3856 #undef __FUNCT__ 3857 #define __FUNCT__ "TSVISetVariableBounds" 3858 /*@ 3859 TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu 3860 3861 Input Parameters: 3862 . ts - the TS context. 3863 . xl - lower bound. 3864 . xu - upper bound. 3865 3866 Notes: 3867 If this routine is not called then the lower and upper bounds are set to 3868 SNES_VI_NINF and SNES_VI_INF respectively during SNESSetUp(). 3869 3870 Level: advanced 3871 3872 @*/ 3873 PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu) 3874 { 3875 PetscErrorCode ierr; 3876 SNES snes; 3877 3878 PetscFunctionBegin; 3879 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3880 ierr = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(ierr); 3881 PetscFunctionReturn(0); 3882 } 3883 3884 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3885 #include <mex.h> 3886 3887 typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext; 3888 3889 #undef __FUNCT__ 3890 #define __FUNCT__ "TSComputeFunction_Matlab" 3891 /* 3892 TSComputeFunction_Matlab - Calls the function that has been set with 3893 TSSetFunctionMatlab(). 3894 3895 Collective on TS 3896 3897 Input Parameters: 3898 + snes - the TS context 3899 - u - input vector 3900 3901 Output Parameter: 3902 . y - function vector, as set by TSSetFunction() 3903 3904 Notes: 3905 TSComputeFunction() is typically used within nonlinear solvers 3906 implementations, so most users would not generally call this routine 3907 themselves. 3908 3909 Level: developer 3910 3911 .keywords: TS, nonlinear, compute, function 3912 3913 .seealso: TSSetFunction(), TSGetFunction() 3914 */ 3915 PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec u,Vec udot,Vec y, void *ctx) 3916 { 3917 PetscErrorCode ierr; 3918 TSMatlabContext *sctx = (TSMatlabContext *)ctx; 3919 int nlhs = 1,nrhs = 7; 3920 mxArray *plhs[1],*prhs[7]; 3921 long long int lx = 0,lxdot = 0,ly = 0,ls = 0; 3922 3923 PetscFunctionBegin; 3924 PetscValidHeaderSpecific(snes,TS_CLASSID,1); 3925 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 3926 PetscValidHeaderSpecific(udot,VEC_CLASSID,4); 3927 PetscValidHeaderSpecific(y,VEC_CLASSID,5); 3928 PetscCheckSameComm(snes,1,u,3); 3929 PetscCheckSameComm(snes,1,y,5); 3930 3931 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 3932 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 3933 ierr = PetscMemcpy(&lxdot,&udot,sizeof(udot));CHKERRQ(ierr); 3934 ierr = PetscMemcpy(&ly,&y,sizeof(u));CHKERRQ(ierr); 3935 prhs[0] = mxCreateDoubleScalar((double)ls); 3936 prhs[1] = mxCreateDoubleScalar(time); 3937 prhs[2] = mxCreateDoubleScalar((double)lx); 3938 prhs[3] = mxCreateDoubleScalar((double)lxdot); 3939 prhs[4] = mxCreateDoubleScalar((double)ly); 3940 prhs[5] = mxCreateString(sctx->funcname); 3941 prhs[6] = sctx->ctx; 3942 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");CHKERRQ(ierr); 3943 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 3944 mxDestroyArray(prhs[0]); 3945 mxDestroyArray(prhs[1]); 3946 mxDestroyArray(prhs[2]); 3947 mxDestroyArray(prhs[3]); 3948 mxDestroyArray(prhs[4]); 3949 mxDestroyArray(prhs[5]); 3950 mxDestroyArray(plhs[0]); 3951 PetscFunctionReturn(0); 3952 } 3953 3954 3955 #undef __FUNCT__ 3956 #define __FUNCT__ "TSSetFunctionMatlab" 3957 /* 3958 TSSetFunctionMatlab - Sets the function evaluation routine and function 3959 vector for use by the TS routines in solving ODEs 3960 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 3961 3962 Logically Collective on TS 3963 3964 Input Parameters: 3965 + ts - the TS context 3966 - func - function evaluation routine 3967 3968 Calling sequence of func: 3969 $ func (TS ts,PetscReal time,Vec u,Vec udot,Vec f,void *ctx); 3970 3971 Level: beginner 3972 3973 .keywords: TS, nonlinear, set, function 3974 3975 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 3976 */ 3977 PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx) 3978 { 3979 PetscErrorCode ierr; 3980 TSMatlabContext *sctx; 3981 3982 PetscFunctionBegin; 3983 /* currently sctx is memory bleed */ 3984 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 3985 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 3986 /* 3987 This should work, but it doesn't 3988 sctx->ctx = ctx; 3989 mexMakeArrayPersistent(sctx->ctx); 3990 */ 3991 sctx->ctx = mxDuplicateArray(ctx); 3992 ierr = TSSetIFunction(ts,PETSC_NULL,TSComputeFunction_Matlab,sctx);CHKERRQ(ierr); 3993 PetscFunctionReturn(0); 3994 } 3995 3996 #undef __FUNCT__ 3997 #define __FUNCT__ "TSComputeJacobian_Matlab" 3998 /* 3999 TSComputeJacobian_Matlab - Calls the function that has been set with 4000 TSSetJacobianMatlab(). 4001 4002 Collective on TS 4003 4004 Input Parameters: 4005 + ts - the TS context 4006 . u - input vector 4007 . A, B - the matrices 4008 - ctx - user context 4009 4010 Output Parameter: 4011 . flag - structure of the matrix 4012 4013 Level: developer 4014 4015 .keywords: TS, nonlinear, compute, function 4016 4017 .seealso: TSSetFunction(), TSGetFunction() 4018 @*/ 4019 PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec u,Vec udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4020 { 4021 PetscErrorCode ierr; 4022 TSMatlabContext *sctx = (TSMatlabContext *)ctx; 4023 int nlhs = 2,nrhs = 9; 4024 mxArray *plhs[2],*prhs[9]; 4025 long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0; 4026 4027 PetscFunctionBegin; 4028 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4029 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 4030 4031 /* call Matlab function in ctx with arguments u and y */ 4032 4033 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 4034 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4035 ierr = PetscMemcpy(&lxdot,&udot,sizeof(u));CHKERRQ(ierr); 4036 ierr = PetscMemcpy(&lA,A,sizeof(u));CHKERRQ(ierr); 4037 ierr = PetscMemcpy(&lB,B,sizeof(u));CHKERRQ(ierr); 4038 prhs[0] = mxCreateDoubleScalar((double)ls); 4039 prhs[1] = mxCreateDoubleScalar((double)time); 4040 prhs[2] = mxCreateDoubleScalar((double)lx); 4041 prhs[3] = mxCreateDoubleScalar((double)lxdot); 4042 prhs[4] = mxCreateDoubleScalar((double)shift); 4043 prhs[5] = mxCreateDoubleScalar((double)lA); 4044 prhs[6] = mxCreateDoubleScalar((double)lB); 4045 prhs[7] = mxCreateString(sctx->funcname); 4046 prhs[8] = sctx->ctx; 4047 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");CHKERRQ(ierr); 4048 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4049 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4050 mxDestroyArray(prhs[0]); 4051 mxDestroyArray(prhs[1]); 4052 mxDestroyArray(prhs[2]); 4053 mxDestroyArray(prhs[3]); 4054 mxDestroyArray(prhs[4]); 4055 mxDestroyArray(prhs[5]); 4056 mxDestroyArray(prhs[6]); 4057 mxDestroyArray(prhs[7]); 4058 mxDestroyArray(plhs[0]); 4059 mxDestroyArray(plhs[1]); 4060 PetscFunctionReturn(0); 4061 } 4062 4063 4064 #undef __FUNCT__ 4065 #define __FUNCT__ "TSSetJacobianMatlab" 4066 /* 4067 TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4068 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 4069 4070 Logically Collective on TS 4071 4072 Input Parameters: 4073 + ts - the TS context 4074 . A,B - Jacobian matrices 4075 . func - function evaluation routine 4076 - ctx - user context 4077 4078 Calling sequence of func: 4079 $ flag = func (TS ts,PetscReal time,Vec u,Vec udot,Mat A,Mat B,void *ctx); 4080 4081 4082 Level: developer 4083 4084 .keywords: TS, nonlinear, set, function 4085 4086 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4087 */ 4088 PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx) 4089 { 4090 PetscErrorCode ierr; 4091 TSMatlabContext *sctx; 4092 4093 PetscFunctionBegin; 4094 /* currently sctx is memory bleed */ 4095 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4096 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4097 /* 4098 This should work, but it doesn't 4099 sctx->ctx = ctx; 4100 mexMakeArrayPersistent(sctx->ctx); 4101 */ 4102 sctx->ctx = mxDuplicateArray(ctx); 4103 ierr = TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4104 PetscFunctionReturn(0); 4105 } 4106 4107 #undef __FUNCT__ 4108 #define __FUNCT__ "TSMonitor_Matlab" 4109 /* 4110 TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab(). 4111 4112 Collective on TS 4113 4114 .seealso: TSSetFunction(), TSGetFunction() 4115 @*/ 4116 PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec u, void *ctx) 4117 { 4118 PetscErrorCode ierr; 4119 TSMatlabContext *sctx = (TSMatlabContext *)ctx; 4120 int nlhs = 1,nrhs = 6; 4121 mxArray *plhs[1],*prhs[6]; 4122 long long int lx = 0,ls = 0; 4123 4124 PetscFunctionBegin; 4125 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4126 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 4127 4128 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 4129 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4130 prhs[0] = mxCreateDoubleScalar((double)ls); 4131 prhs[1] = mxCreateDoubleScalar((double)it); 4132 prhs[2] = mxCreateDoubleScalar((double)time); 4133 prhs[3] = mxCreateDoubleScalar((double)lx); 4134 prhs[4] = mxCreateString(sctx->funcname); 4135 prhs[5] = sctx->ctx; 4136 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");CHKERRQ(ierr); 4137 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4138 mxDestroyArray(prhs[0]); 4139 mxDestroyArray(prhs[1]); 4140 mxDestroyArray(prhs[2]); 4141 mxDestroyArray(prhs[3]); 4142 mxDestroyArray(prhs[4]); 4143 mxDestroyArray(plhs[0]); 4144 PetscFunctionReturn(0); 4145 } 4146 4147 4148 #undef __FUNCT__ 4149 #define __FUNCT__ "TSMonitorSetMatlab" 4150 /* 4151 TSMonitorSetMatlab - Sets the monitor function from Matlab 4152 4153 Level: developer 4154 4155 .keywords: TS, nonlinear, set, function 4156 4157 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4158 */ 4159 PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx) 4160 { 4161 PetscErrorCode ierr; 4162 TSMatlabContext *sctx; 4163 4164 PetscFunctionBegin; 4165 /* currently sctx is memory bleed */ 4166 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4167 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4168 /* 4169 This should work, but it doesn't 4170 sctx->ctx = ctx; 4171 mexMakeArrayPersistent(sctx->ctx); 4172 */ 4173 sctx->ctx = mxDuplicateArray(ctx); 4174 ierr = TSMonitorSet(ts,TSMonitor_Matlab,sctx,PETSC_NULL);CHKERRQ(ierr); 4175 PetscFunctionReturn(0); 4176 } 4177 #endif 4178 4179 4180 4181 #undef __FUNCT__ 4182 #define __FUNCT__ "TSMonitorLGSolution" 4183 /*@C 4184 TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector 4185 in a time based line graph 4186 4187 Collective on TS 4188 4189 Input Parameters: 4190 + ts - the TS context 4191 . step - current time-step 4192 . ptime - current time 4193 - lg - a line graph object 4194 4195 Level: intermediate 4196 4197 Notes: each process in a parallel run displays its component solutions in a separate window 4198 4199 .keywords: TS, vector, monitor, view 4200 4201 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4202 @*/ 4203 PetscErrorCode TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4204 { 4205 PetscErrorCode ierr; 4206 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 4207 const PetscScalar *yy; 4208 PetscInt dim; 4209 4210 PetscFunctionBegin; 4211 if (!step) { 4212 PetscDrawAxis axis; 4213 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4214 ierr = PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");CHKERRQ(ierr); 4215 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 4216 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 4217 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4218 } 4219 ierr = VecGetArrayRead(u,&yy);CHKERRQ(ierr); 4220 #if defined(PETSC_USE_COMPLEX) 4221 { 4222 PetscReal *yreal; 4223 PetscInt i,n; 4224 ierr = VecGetLocalSize(u,&n);CHKERRQ(ierr); 4225 ierr = PetscMalloc(n*sizeof(PetscReal),&yreal);CHKERRQ(ierr); 4226 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 4227 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 4228 ierr = PetscFree(yreal);CHKERRQ(ierr); 4229 } 4230 #else 4231 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 4232 #endif 4233 ierr = VecRestoreArrayRead(u,&yy);CHKERRQ(ierr); 4234 if (((ctx->howoften > 0) && (!(step % ctx->howoften)) && (step > -1)) || ((ctx->howoften == -1) && (step == -1))){ 4235 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4236 } 4237 PetscFunctionReturn(0); 4238 } 4239 4240 #undef __FUNCT__ 4241 #define __FUNCT__ "TSMonitorLGError" 4242 /*@C 4243 TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the solution vector 4244 in a time based line graph 4245 4246 Collective on TS 4247 4248 Input Parameters: 4249 + ts - the TS context 4250 . step - current time-step 4251 . ptime - current time 4252 - lg - a line graph object 4253 4254 Level: intermediate 4255 4256 Notes: 4257 Only for sequential solves. 4258 4259 The user must provide the solution using TSSetSolutionFunction() to use this monitor. 4260 4261 Options Database Keys: 4262 . -ts_monitor_lg_error - create a graphical monitor of error history 4263 4264 .keywords: TS, vector, monitor, view 4265 4266 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction() 4267 @*/ 4268 PetscErrorCode TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4269 { 4270 PetscErrorCode ierr; 4271 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 4272 const PetscScalar *yy; 4273 Vec y; 4274 PetscInt dim; 4275 4276 PetscFunctionBegin; 4277 if (!step) { 4278 PetscDrawAxis axis; 4279 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4280 ierr = PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Solution");CHKERRQ(ierr); 4281 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 4282 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 4283 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4284 } 4285 ierr = VecDuplicate(u,&y);CHKERRQ(ierr); 4286 ierr = TSComputeSolutionFunction(ts,ptime,y);CHKERRQ(ierr); 4287 ierr = VecAXPY(y,-1.0,u);CHKERRQ(ierr); 4288 ierr = VecGetArrayRead(y,&yy);CHKERRQ(ierr); 4289 #if defined(PETSC_USE_COMPLEX) 4290 { 4291 PetscReal *yreal; 4292 PetscInt i,n; 4293 ierr = VecGetLocalSize(y,&n);CHKERRQ(ierr); 4294 ierr = PetscMalloc(n*sizeof(PetscReal),&yreal);CHKERRQ(ierr); 4295 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 4296 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 4297 ierr = PetscFree(yreal);CHKERRQ(ierr); 4298 } 4299 #else 4300 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 4301 #endif 4302 ierr = VecRestoreArrayRead(y,&yy);CHKERRQ(ierr); 4303 ierr = VecDestroy(&y);CHKERRQ(ierr); 4304 if (((ctx->howoften > 0) && (!(step % ctx->howoften)) && (step > -1)) || ((ctx->howoften == -1) && (step == -1))){ 4305 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4306 } 4307 PetscFunctionReturn(0); 4308 } 4309 4310 #undef __FUNCT__ 4311 #define __FUNCT__ "TSMonitorLGSNESIterations" 4312 PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 4313 { 4314 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 4315 PetscReal x = ptime,y; 4316 PetscErrorCode ierr; 4317 PetscInt its; 4318 4319 PetscFunctionBegin; 4320 if (!n) { 4321 PetscDrawAxis axis; 4322 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4323 ierr = PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");CHKERRQ(ierr); 4324 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4325 ctx->snes_its = 0; 4326 } 4327 ierr = TSGetSNESIterations(ts,&its);CHKERRQ(ierr); 4328 y = its - ctx->snes_its; 4329 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 4330 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))){ 4331 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4332 } 4333 ctx->snes_its = its; 4334 PetscFunctionReturn(0); 4335 } 4336 4337 #undef __FUNCT__ 4338 #define __FUNCT__ "TSMonitorLGKSPIterations" 4339 PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 4340 { 4341 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 4342 PetscReal x = ptime,y; 4343 PetscErrorCode ierr; 4344 PetscInt its; 4345 4346 PetscFunctionBegin; 4347 if (!n) { 4348 PetscDrawAxis axis; 4349 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4350 ierr = PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");CHKERRQ(ierr); 4351 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4352 ctx->ksp_its = 0; 4353 } 4354 ierr = TSGetKSPIterations(ts,&its);CHKERRQ(ierr); 4355 y = its - ctx->ksp_its; 4356 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 4357 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))){ 4358 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4359 } 4360 ctx->ksp_its = its; 4361 PetscFunctionReturn(0); 4362 } 4363 4364 #undef __FUNCT__ 4365 #define __FUNCT__ "TSComputeLinearStability" 4366 /*@ 4367 TSComputeLinearStability - computes the linear stability function at a point 4368 4369 Collective on TS and Vec 4370 4371 Input Parameters: 4372 + ts - the TS context 4373 - xr,xi - real and imaginary part of input arguments 4374 4375 Output Parameters: 4376 . yr,yi - real and imaginary part of function value 4377 4378 Level: developer 4379 4380 .keywords: TS, compute 4381 4382 .seealso: TSSetRHSFunction(), TSComputeIFunction() 4383 @*/ 4384 PetscErrorCode TSComputeLinearStability(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi) 4385 { 4386 PetscErrorCode ierr; 4387 4388 PetscFunctionBegin; 4389 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4390 if (!ts->ops->linearstability) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_SUP,"Linearized stability function not provided for this method"); 4391 ierr = (*ts->ops->linearstability)(ts,xr,xi,yr,yi);CHKERRQ(ierr); 4392 PetscFunctionReturn(0); 4393 } 4394