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