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