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