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 /* ----- Routines to initialize and destroy a timestepper ---- */ 1731 #undef __FUNCT__ 1732 #define __FUNCT__ "TSSetProblemType" 1733 /*@ 1734 TSSetProblemType - Sets the type of problem to be solved. 1735 1736 Not collective 1737 1738 Input Parameters: 1739 + ts - The TS 1740 - type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms 1741 .vb 1742 U_t - A U = 0 (linear) 1743 U_t - A(t) U = 0 (linear) 1744 F(t,U,U_t) = 0 (nonlinear) 1745 .ve 1746 1747 Level: beginner 1748 1749 .keywords: TS, problem type 1750 .seealso: TSSetUp(), TSProblemType, TS 1751 @*/ 1752 PetscErrorCode TSSetProblemType(TS ts, TSProblemType type) 1753 { 1754 PetscErrorCode ierr; 1755 1756 PetscFunctionBegin; 1757 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1758 ts->problem_type = type; 1759 if (type == TS_LINEAR) { 1760 SNES snes; 1761 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1762 ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 1763 } 1764 PetscFunctionReturn(0); 1765 } 1766 1767 #undef __FUNCT__ 1768 #define __FUNCT__ "TSGetProblemType" 1769 /*@C 1770 TSGetProblemType - Gets the type of problem to be solved. 1771 1772 Not collective 1773 1774 Input Parameter: 1775 . ts - The TS 1776 1777 Output Parameter: 1778 . type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms 1779 .vb 1780 M U_t = A U 1781 M(t) U_t = A(t) U 1782 F(t,U,U_t) 1783 .ve 1784 1785 Level: beginner 1786 1787 .keywords: TS, problem type 1788 .seealso: TSSetUp(), TSProblemType, TS 1789 @*/ 1790 PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type) 1791 { 1792 PetscFunctionBegin; 1793 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1794 PetscValidIntPointer(type,2); 1795 *type = ts->problem_type; 1796 PetscFunctionReturn(0); 1797 } 1798 1799 #undef __FUNCT__ 1800 #define __FUNCT__ "TSSetUp" 1801 /*@ 1802 TSSetUp - Sets up the internal data structures for the later use 1803 of a timestepper. 1804 1805 Collective on TS 1806 1807 Input Parameter: 1808 . ts - the TS context obtained from TSCreate() 1809 1810 Notes: 1811 For basic use of the TS solvers the user need not explicitly call 1812 TSSetUp(), since these actions will automatically occur during 1813 the call to TSStep(). However, if one wishes to control this 1814 phase separately, TSSetUp() should be called after TSCreate() 1815 and optional routines of the form TSSetXXX(), but before TSStep(). 1816 1817 Level: advanced 1818 1819 .keywords: TS, timestep, setup 1820 1821 .seealso: TSCreate(), TSStep(), TSDestroy() 1822 @*/ 1823 PetscErrorCode TSSetUp(TS ts) 1824 { 1825 PetscErrorCode ierr; 1826 DM dm; 1827 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 1828 PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*); 1829 TSIJacobian ijac; 1830 TSRHSJacobian rhsjac; 1831 1832 PetscFunctionBegin; 1833 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1834 if (ts->setupcalled) PetscFunctionReturn(0); 1835 1836 if (!((PetscObject)ts)->type_name) { 1837 ierr = TSSetType(ts,TSEULER);CHKERRQ(ierr); 1838 } 1839 1840 if (!ts->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetSolution() first"); 1841 1842 ierr = TSGetAdapt(ts,&ts->adapt);CHKERRQ(ierr); 1843 1844 if (ts->rhsjacobian.reuse) { 1845 Mat Amat,Pmat; 1846 SNES snes; 1847 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1848 ierr = SNESGetJacobian(snes,&Amat,&Pmat,NULL,NULL);CHKERRQ(ierr); 1849 /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would 1850 * have displaced the RHS matrix */ 1851 if (Amat == ts->Arhs) { 1852 ierr = MatDuplicate(ts->Arhs,MAT_DO_NOT_COPY_VALUES,&Amat);CHKERRQ(ierr); 1853 ierr = SNESSetJacobian(snes,Amat,NULL,NULL,NULL);CHKERRQ(ierr); 1854 ierr = MatDestroy(&Amat);CHKERRQ(ierr); 1855 } 1856 if (Pmat == ts->Brhs) { 1857 ierr = MatDuplicate(ts->Brhs,MAT_DO_NOT_COPY_VALUES,&Pmat);CHKERRQ(ierr); 1858 ierr = SNESSetJacobian(snes,NULL,Pmat,NULL,NULL);CHKERRQ(ierr); 1859 ierr = MatDestroy(&Pmat);CHKERRQ(ierr); 1860 } 1861 } 1862 1863 if (ts->ops->setup) { 1864 ierr = (*ts->ops->setup)(ts);CHKERRQ(ierr); 1865 } 1866 1867 /* in the case where we've set a DMTSFunction or what have you, we need the default SNESFunction 1868 to be set right but can't do it elsewhere due to the overreliance on ctx=ts. 1869 */ 1870 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1871 ierr = DMSNESGetFunction(dm,&func,NULL);CHKERRQ(ierr); 1872 if (!func) { 1873 ierr =DMSNESSetFunction(dm,SNESTSFormFunction,ts);CHKERRQ(ierr); 1874 } 1875 /* if the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it. 1876 Otherwise, the SNES will use coloring internally to form the Jacobian. 1877 */ 1878 ierr = DMSNESGetJacobian(dm,&jac,NULL);CHKERRQ(ierr); 1879 ierr = DMTSGetIJacobian(dm,&ijac,NULL);CHKERRQ(ierr); 1880 ierr = DMTSGetRHSJacobian(dm,&rhsjac,NULL);CHKERRQ(ierr); 1881 if (!jac && (ijac || rhsjac)) { 1882 ierr = DMSNESSetJacobian(dm,SNESTSFormJacobian,ts);CHKERRQ(ierr); 1883 } 1884 ts->setupcalled = PETSC_TRUE; 1885 PetscFunctionReturn(0); 1886 } 1887 1888 #undef __FUNCT__ 1889 #define __FUNCT__ "TSReset" 1890 /*@ 1891 TSReset - Resets a TS context and removes any allocated Vecs and Mats. 1892 1893 Collective on TS 1894 1895 Input Parameter: 1896 . ts - the TS context obtained from TSCreate() 1897 1898 Level: beginner 1899 1900 .keywords: TS, timestep, reset 1901 1902 .seealso: TSCreate(), TSSetup(), TSDestroy() 1903 @*/ 1904 PetscErrorCode TSReset(TS ts) 1905 { 1906 PetscErrorCode ierr; 1907 1908 PetscFunctionBegin; 1909 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1910 if (ts->ops->reset) { 1911 ierr = (*ts->ops->reset)(ts);CHKERRQ(ierr); 1912 } 1913 if (ts->snes) {ierr = SNESReset(ts->snes);CHKERRQ(ierr);} 1914 1915 ierr = MatDestroy(&ts->Arhs);CHKERRQ(ierr); 1916 ierr = MatDestroy(&ts->Brhs);CHKERRQ(ierr); 1917 ierr = VecDestroy(&ts->Frhs);CHKERRQ(ierr); 1918 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 1919 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 1920 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 1921 ierr = VecDestroyVecs(ts->nwork,&ts->work);CHKERRQ(ierr); 1922 1923 ts->setupcalled = PETSC_FALSE; 1924 PetscFunctionReturn(0); 1925 } 1926 1927 #undef __FUNCT__ 1928 #define __FUNCT__ "TSDestroy" 1929 /*@ 1930 TSDestroy - Destroys the timestepper context that was created 1931 with TSCreate(). 1932 1933 Collective on TS 1934 1935 Input Parameter: 1936 . ts - the TS context obtained from TSCreate() 1937 1938 Level: beginner 1939 1940 .keywords: TS, timestepper, destroy 1941 1942 .seealso: TSCreate(), TSSetUp(), TSSolve() 1943 @*/ 1944 PetscErrorCode TSDestroy(TS *ts) 1945 { 1946 PetscErrorCode ierr; 1947 1948 PetscFunctionBegin; 1949 if (!*ts) PetscFunctionReturn(0); 1950 PetscValidHeaderSpecific((*ts),TS_CLASSID,1); 1951 if (--((PetscObject)(*ts))->refct > 0) {*ts = 0; PetscFunctionReturn(0);} 1952 1953 ierr = TSReset((*ts));CHKERRQ(ierr); 1954 1955 /* if memory was published with SAWs then destroy it */ 1956 ierr = PetscObjectSAWsViewOff((PetscObject)*ts);CHKERRQ(ierr); 1957 if ((*ts)->ops->destroy) {ierr = (*(*ts)->ops->destroy)((*ts));CHKERRQ(ierr);} 1958 1959 ierr = TSAdaptDestroy(&(*ts)->adapt);CHKERRQ(ierr); 1960 if ((*ts)->event) { 1961 ierr = TSEventMonitorDestroy(&(*ts)->event);CHKERRQ(ierr); 1962 } 1963 ierr = SNESDestroy(&(*ts)->snes);CHKERRQ(ierr); 1964 ierr = DMDestroy(&(*ts)->dm);CHKERRQ(ierr); 1965 ierr = TSMonitorCancel((*ts));CHKERRQ(ierr); 1966 1967 ierr = PetscHeaderDestroy(ts);CHKERRQ(ierr); 1968 PetscFunctionReturn(0); 1969 } 1970 1971 #undef __FUNCT__ 1972 #define __FUNCT__ "TSGetSNES" 1973 /*@ 1974 TSGetSNES - Returns the SNES (nonlinear solver) associated with 1975 a TS (timestepper) context. Valid only for nonlinear problems. 1976 1977 Not Collective, but SNES is parallel if TS is parallel 1978 1979 Input Parameter: 1980 . ts - the TS context obtained from TSCreate() 1981 1982 Output Parameter: 1983 . snes - the nonlinear solver context 1984 1985 Notes: 1986 The user can then directly manipulate the SNES context to set various 1987 options, etc. Likewise, the user can then extract and manipulate the 1988 KSP, KSP, and PC contexts as well. 1989 1990 TSGetSNES() does not work for integrators that do not use SNES; in 1991 this case TSGetSNES() returns NULL in snes. 1992 1993 Level: beginner 1994 1995 .keywords: timestep, get, SNES 1996 @*/ 1997 PetscErrorCode TSGetSNES(TS ts,SNES *snes) 1998 { 1999 PetscErrorCode ierr; 2000 2001 PetscFunctionBegin; 2002 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2003 PetscValidPointer(snes,2); 2004 if (!ts->snes) { 2005 ierr = SNESCreate(PetscObjectComm((PetscObject)ts),&ts->snes);CHKERRQ(ierr); 2006 ierr = SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 2007 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->snes);CHKERRQ(ierr); 2008 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);CHKERRQ(ierr); 2009 if (ts->dm) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 2010 if (ts->problem_type == TS_LINEAR) { 2011 ierr = SNESSetType(ts->snes,SNESKSPONLY);CHKERRQ(ierr); 2012 } 2013 } 2014 *snes = ts->snes; 2015 PetscFunctionReturn(0); 2016 } 2017 2018 #undef __FUNCT__ 2019 #define __FUNCT__ "TSSetSNES" 2020 /*@ 2021 TSSetSNES - Set the SNES (nonlinear solver) to be used by the timestepping context 2022 2023 Collective 2024 2025 Input Parameter: 2026 + ts - the TS context obtained from TSCreate() 2027 - snes - the nonlinear solver context 2028 2029 Notes: 2030 Most users should have the TS created by calling TSGetSNES() 2031 2032 Level: developer 2033 2034 .keywords: timestep, set, SNES 2035 @*/ 2036 PetscErrorCode TSSetSNES(TS ts,SNES snes) 2037 { 2038 PetscErrorCode ierr; 2039 PetscErrorCode (*func)(SNES,Vec,Mat,Mat,void*); 2040 2041 PetscFunctionBegin; 2042 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2043 PetscValidHeaderSpecific(snes,SNES_CLASSID,2); 2044 ierr = PetscObjectReference((PetscObject)snes);CHKERRQ(ierr); 2045 ierr = SNESDestroy(&ts->snes);CHKERRQ(ierr); 2046 2047 ts->snes = snes; 2048 2049 ierr = SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 2050 ierr = SNESGetJacobian(ts->snes,NULL,NULL,&func,NULL);CHKERRQ(ierr); 2051 if (func == SNESTSFormJacobian) { 2052 ierr = SNESSetJacobian(ts->snes,NULL,NULL,SNESTSFormJacobian,ts);CHKERRQ(ierr); 2053 } 2054 PetscFunctionReturn(0); 2055 } 2056 2057 #undef __FUNCT__ 2058 #define __FUNCT__ "TSGetKSP" 2059 /*@ 2060 TSGetKSP - Returns the KSP (linear solver) associated with 2061 a TS (timestepper) context. 2062 2063 Not Collective, but KSP is parallel if TS is parallel 2064 2065 Input Parameter: 2066 . ts - the TS context obtained from TSCreate() 2067 2068 Output Parameter: 2069 . ksp - the nonlinear solver context 2070 2071 Notes: 2072 The user can then directly manipulate the KSP context to set various 2073 options, etc. Likewise, the user can then extract and manipulate the 2074 KSP and PC contexts as well. 2075 2076 TSGetKSP() does not work for integrators that do not use KSP; 2077 in this case TSGetKSP() returns NULL in ksp. 2078 2079 Level: beginner 2080 2081 .keywords: timestep, get, KSP 2082 @*/ 2083 PetscErrorCode TSGetKSP(TS ts,KSP *ksp) 2084 { 2085 PetscErrorCode ierr; 2086 SNES snes; 2087 2088 PetscFunctionBegin; 2089 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2090 PetscValidPointer(ksp,2); 2091 if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first"); 2092 if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()"); 2093 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2094 ierr = SNESGetKSP(snes,ksp);CHKERRQ(ierr); 2095 PetscFunctionReturn(0); 2096 } 2097 2098 /* ----------- Routines to set solver parameters ---------- */ 2099 2100 #undef __FUNCT__ 2101 #define __FUNCT__ "TSGetDuration" 2102 /*@ 2103 TSGetDuration - Gets the maximum number of timesteps to use and 2104 maximum time for iteration. 2105 2106 Not Collective 2107 2108 Input Parameters: 2109 + ts - the TS context obtained from TSCreate() 2110 . maxsteps - maximum number of iterations to use, or NULL 2111 - maxtime - final time to iterate to, or NULL 2112 2113 Level: intermediate 2114 2115 .keywords: TS, timestep, get, maximum, iterations, time 2116 @*/ 2117 PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime) 2118 { 2119 PetscFunctionBegin; 2120 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2121 if (maxsteps) { 2122 PetscValidIntPointer(maxsteps,2); 2123 *maxsteps = ts->max_steps; 2124 } 2125 if (maxtime) { 2126 PetscValidScalarPointer(maxtime,3); 2127 *maxtime = ts->max_time; 2128 } 2129 PetscFunctionReturn(0); 2130 } 2131 2132 #undef __FUNCT__ 2133 #define __FUNCT__ "TSSetDuration" 2134 /*@ 2135 TSSetDuration - Sets the maximum number of timesteps to use and 2136 maximum time for iteration. 2137 2138 Logically Collective on TS 2139 2140 Input Parameters: 2141 + ts - the TS context obtained from TSCreate() 2142 . maxsteps - maximum number of iterations to use 2143 - maxtime - final time to iterate to 2144 2145 Options Database Keys: 2146 . -ts_max_steps <maxsteps> - Sets maxsteps 2147 . -ts_final_time <maxtime> - Sets maxtime 2148 2149 Notes: 2150 The default maximum number of iterations is 5000. Default time is 5.0 2151 2152 Level: intermediate 2153 2154 .keywords: TS, timestep, set, maximum, iterations 2155 2156 .seealso: TSSetExactFinalTime() 2157 @*/ 2158 PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime) 2159 { 2160 PetscFunctionBegin; 2161 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2162 PetscValidLogicalCollectiveInt(ts,maxsteps,2); 2163 PetscValidLogicalCollectiveReal(ts,maxtime,2); 2164 if (maxsteps >= 0) ts->max_steps = maxsteps; 2165 if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime; 2166 PetscFunctionReturn(0); 2167 } 2168 2169 #undef __FUNCT__ 2170 #define __FUNCT__ "TSSetSolution" 2171 /*@ 2172 TSSetSolution - Sets the initial solution vector 2173 for use by the TS routines. 2174 2175 Logically Collective on TS and Vec 2176 2177 Input Parameters: 2178 + ts - the TS context obtained from TSCreate() 2179 - u - the solution vector 2180 2181 Level: beginner 2182 2183 .keywords: TS, timestep, set, solution, initial conditions 2184 @*/ 2185 PetscErrorCode TSSetSolution(TS ts,Vec u) 2186 { 2187 PetscErrorCode ierr; 2188 DM dm; 2189 2190 PetscFunctionBegin; 2191 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2192 PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2193 ierr = PetscObjectReference((PetscObject)u);CHKERRQ(ierr); 2194 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 2195 2196 ts->vec_sol = u; 2197 2198 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 2199 ierr = DMShellSetGlobalVector(dm,u);CHKERRQ(ierr); 2200 PetscFunctionReturn(0); 2201 } 2202 2203 #undef __FUNCT__ 2204 #define __FUNCT__ "TSSetPreStep" 2205 /*@C 2206 TSSetPreStep - Sets the general-purpose function 2207 called once at the beginning of each time step. 2208 2209 Logically Collective on TS 2210 2211 Input Parameters: 2212 + ts - The TS context obtained from TSCreate() 2213 - func - The function 2214 2215 Calling sequence of func: 2216 . func (TS ts); 2217 2218 Level: intermediate 2219 2220 Note: 2221 If a step is rejected, TSStep() will call this routine again before each attempt. 2222 The last completed time step number can be queried using TSGetTimeStepNumber(), the 2223 size of the step being attempted can be obtained using TSGetTimeStep(). 2224 2225 .keywords: TS, timestep 2226 .seealso: TSSetPreStage(), TSSetPostStage(), TSSetPostStep(), TSStep() 2227 @*/ 2228 PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS)) 2229 { 2230 PetscFunctionBegin; 2231 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2232 ts->prestep = func; 2233 PetscFunctionReturn(0); 2234 } 2235 2236 #undef __FUNCT__ 2237 #define __FUNCT__ "TSPreStep" 2238 /*@ 2239 TSPreStep - Runs the user-defined pre-step function. 2240 2241 Collective on TS 2242 2243 Input Parameters: 2244 . ts - The TS context obtained from TSCreate() 2245 2246 Notes: 2247 TSPreStep() is typically used within time stepping implementations, 2248 so most users would not generally call this routine themselves. 2249 2250 Level: developer 2251 2252 .keywords: TS, timestep 2253 .seealso: TSSetPreStep(), TSPreStage(), TSPostStage(), TSPostStep() 2254 @*/ 2255 PetscErrorCode TSPreStep(TS ts) 2256 { 2257 PetscErrorCode ierr; 2258 2259 PetscFunctionBegin; 2260 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2261 if (ts->prestep) { 2262 PetscStackCallStandard((*ts->prestep),(ts)); 2263 } 2264 PetscFunctionReturn(0); 2265 } 2266 2267 #undef __FUNCT__ 2268 #define __FUNCT__ "TSSetPreStage" 2269 /*@C 2270 TSSetPreStage - Sets the general-purpose function 2271 called once at the beginning of each stage. 2272 2273 Logically Collective on TS 2274 2275 Input Parameters: 2276 + ts - The TS context obtained from TSCreate() 2277 - func - The function 2278 2279 Calling sequence of func: 2280 . PetscErrorCode func(TS ts, PetscReal stagetime); 2281 2282 Level: intermediate 2283 2284 Note: 2285 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 2286 The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being 2287 attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime(). 2288 2289 .keywords: TS, timestep 2290 .seealso: TSSetPostStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 2291 @*/ 2292 PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS,PetscReal)) 2293 { 2294 PetscFunctionBegin; 2295 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2296 ts->prestage = func; 2297 PetscFunctionReturn(0); 2298 } 2299 2300 #undef __FUNCT__ 2301 #define __FUNCT__ "TSSetPostStage" 2302 /*@C 2303 TSSetPostStage - Sets the general-purpose function 2304 called once at the end of each stage. 2305 2306 Logically Collective on TS 2307 2308 Input Parameters: 2309 + ts - The TS context obtained from TSCreate() 2310 - func - The function 2311 2312 Calling sequence of func: 2313 . PetscErrorCode func(TS ts, PetscReal stagetime, PetscInt stageindex, Vec* Y); 2314 2315 Level: intermediate 2316 2317 Note: 2318 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 2319 The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being 2320 attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime(). 2321 2322 .keywords: TS, timestep 2323 .seealso: TSSetPreStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 2324 @*/ 2325 PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS,PetscReal,PetscInt,Vec*)) 2326 { 2327 PetscFunctionBegin; 2328 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2329 ts->poststage = func; 2330 PetscFunctionReturn(0); 2331 } 2332 2333 #undef __FUNCT__ 2334 #define __FUNCT__ "TSPreStage" 2335 /*@ 2336 TSPreStage - Runs the user-defined pre-stage function set using TSSetPreStage() 2337 2338 Collective on TS 2339 2340 Input Parameters: 2341 . ts - The TS context obtained from TSCreate() 2342 stagetime - The absolute time of the current stage 2343 2344 Notes: 2345 TSPreStage() is typically used within time stepping implementations, 2346 most users would not generally call this routine themselves. 2347 2348 Level: developer 2349 2350 .keywords: TS, timestep 2351 .seealso: TSPostStage(), TSSetPreStep(), TSPreStep(), TSPostStep() 2352 @*/ 2353 PetscErrorCode TSPreStage(TS ts, PetscReal stagetime) 2354 { 2355 PetscErrorCode ierr; 2356 2357 PetscFunctionBegin; 2358 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2359 if (ts->prestage) { 2360 PetscStackCallStandard((*ts->prestage),(ts,stagetime)); 2361 } 2362 PetscFunctionReturn(0); 2363 } 2364 2365 #undef __FUNCT__ 2366 #define __FUNCT__ "TSPostStage" 2367 /*@ 2368 TSPostStage - Runs the user-defined post-stage function set using TSSetPostStage() 2369 2370 Collective on TS 2371 2372 Input Parameters: 2373 . ts - The TS context obtained from TSCreate() 2374 stagetime - The absolute time of the current stage 2375 stageindex - Stage number 2376 Y - Array of vectors (of size = total number 2377 of stages) with the stage solutions 2378 2379 Notes: 2380 TSPostStage() is typically used within time stepping implementations, 2381 most users would not generally call this routine themselves. 2382 2383 Level: developer 2384 2385 .keywords: TS, timestep 2386 .seealso: TSPreStage(), TSSetPreStep(), TSPreStep(), TSPostStep() 2387 @*/ 2388 PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y) 2389 { 2390 PetscErrorCode ierr; 2391 2392 PetscFunctionBegin; 2393 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2394 if (ts->poststage) { 2395 PetscStackCallStandard((*ts->poststage),(ts,stagetime,stageindex,Y)); 2396 } 2397 PetscFunctionReturn(0); 2398 } 2399 2400 #undef __FUNCT__ 2401 #define __FUNCT__ "TSSetPostStep" 2402 /*@C 2403 TSSetPostStep - Sets the general-purpose function 2404 called once at the end of each time step. 2405 2406 Logically Collective on TS 2407 2408 Input Parameters: 2409 + ts - The TS context obtained from TSCreate() 2410 - func - The function 2411 2412 Calling sequence of func: 2413 $ func (TS ts); 2414 2415 Level: intermediate 2416 2417 .keywords: TS, timestep 2418 .seealso: TSSetPreStep(), TSSetPreStage(), TSGetTimeStep(), TSGetTimeStepNumber(), TSGetTime() 2419 @*/ 2420 PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS)) 2421 { 2422 PetscFunctionBegin; 2423 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2424 ts->poststep = func; 2425 PetscFunctionReturn(0); 2426 } 2427 2428 #undef __FUNCT__ 2429 #define __FUNCT__ "TSPostStep" 2430 /*@ 2431 TSPostStep - Runs the user-defined post-step function. 2432 2433 Collective on TS 2434 2435 Input Parameters: 2436 . ts - The TS context obtained from TSCreate() 2437 2438 Notes: 2439 TSPostStep() is typically used within time stepping implementations, 2440 so most users would not generally call this routine themselves. 2441 2442 Level: developer 2443 2444 .keywords: TS, timestep 2445 @*/ 2446 PetscErrorCode TSPostStep(TS ts) 2447 { 2448 PetscErrorCode ierr; 2449 2450 PetscFunctionBegin; 2451 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2452 if (ts->poststep) { 2453 PetscStackCallStandard((*ts->poststep),(ts)); 2454 } 2455 PetscFunctionReturn(0); 2456 } 2457 2458 /* ------------ Routines to set performance monitoring options ----------- */ 2459 2460 #undef __FUNCT__ 2461 #define __FUNCT__ "TSMonitorSet" 2462 /*@C 2463 TSMonitorSet - Sets an ADDITIONAL function that is to be used at every 2464 timestep to display the iteration's progress. 2465 2466 Logically Collective on TS 2467 2468 Input Parameters: 2469 + ts - the TS context obtained from TSCreate() 2470 . monitor - monitoring routine 2471 . mctx - [optional] user-defined context for private data for the 2472 monitor routine (use NULL if no context is desired) 2473 - monitordestroy - [optional] routine that frees monitor context 2474 (may be NULL) 2475 2476 Calling sequence of monitor: 2477 $ int monitor(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx) 2478 2479 + ts - the TS context 2480 . 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 2481 been interpolated to) 2482 . time - current time 2483 . u - current iterate 2484 - mctx - [optional] monitoring context 2485 2486 Notes: 2487 This routine adds an additional monitor to the list of monitors that 2488 already has been loaded. 2489 2490 Fortran notes: Only a single monitor function can be set for each TS object 2491 2492 Level: intermediate 2493 2494 .keywords: TS, timestep, set, monitor 2495 2496 .seealso: TSMonitorDefault(), TSMonitorCancel() 2497 @*/ 2498 PetscErrorCode TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**)) 2499 { 2500 PetscFunctionBegin; 2501 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2502 if (ts->numbermonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 2503 ts->monitor[ts->numbermonitors] = monitor; 2504 ts->monitordestroy[ts->numbermonitors] = mdestroy; 2505 ts->monitorcontext[ts->numbermonitors++] = (void*)mctx; 2506 PetscFunctionReturn(0); 2507 } 2508 2509 #undef __FUNCT__ 2510 #define __FUNCT__ "TSMonitorCancel" 2511 /*@C 2512 TSMonitorCancel - Clears all the monitors that have been set on a time-step object. 2513 2514 Logically Collective on TS 2515 2516 Input Parameters: 2517 . ts - the TS context obtained from TSCreate() 2518 2519 Notes: 2520 There is no way to remove a single, specific monitor. 2521 2522 Level: intermediate 2523 2524 .keywords: TS, timestep, set, monitor 2525 2526 .seealso: TSMonitorDefault(), TSMonitorSet() 2527 @*/ 2528 PetscErrorCode TSMonitorCancel(TS ts) 2529 { 2530 PetscErrorCode ierr; 2531 PetscInt i; 2532 2533 PetscFunctionBegin; 2534 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2535 for (i=0; i<ts->numbermonitors; i++) { 2536 if (ts->monitordestroy[i]) { 2537 ierr = (*ts->monitordestroy[i])(&ts->monitorcontext[i]);CHKERRQ(ierr); 2538 } 2539 } 2540 ts->numbermonitors = 0; 2541 PetscFunctionReturn(0); 2542 } 2543 2544 #undef __FUNCT__ 2545 #define __FUNCT__ "TSMonitorDefault" 2546 /*@ 2547 TSMonitorDefault - Sets the Default monitor 2548 2549 Level: intermediate 2550 2551 .keywords: TS, set, monitor 2552 2553 .seealso: TSMonitorDefault(), TSMonitorSet() 2554 @*/ 2555 PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy) 2556 { 2557 PetscErrorCode ierr; 2558 PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ts)); 2559 2560 PetscFunctionBegin; 2561 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 2562 ierr = PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g\n",step,(double)ts->time_step,(double)ptime);CHKERRQ(ierr); 2563 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 2564 PetscFunctionReturn(0); 2565 } 2566 2567 #undef __FUNCT__ 2568 #define __FUNCT__ "TSSetRetainStages" 2569 /*@ 2570 TSSetRetainStages - Request that all stages in the upcoming step be stored so that interpolation will be available. 2571 2572 Logically Collective on TS 2573 2574 Input Argument: 2575 . ts - time stepping context 2576 2577 Output Argument: 2578 . flg - PETSC_TRUE or PETSC_FALSE 2579 2580 Level: intermediate 2581 2582 .keywords: TS, set 2583 2584 .seealso: TSInterpolate(), TSSetPostStep() 2585 @*/ 2586 PetscErrorCode TSSetRetainStages(TS ts,PetscBool flg) 2587 { 2588 PetscFunctionBegin; 2589 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2590 ts->retain_stages = flg; 2591 PetscFunctionReturn(0); 2592 } 2593 2594 #undef __FUNCT__ 2595 #define __FUNCT__ "TSInterpolate" 2596 /*@ 2597 TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval 2598 2599 Collective on TS 2600 2601 Input Argument: 2602 + ts - time stepping context 2603 - t - time to interpolate to 2604 2605 Output Argument: 2606 . U - state at given time 2607 2608 Notes: 2609 The user should call TSSetRetainStages() before taking a step in which interpolation will be requested. 2610 2611 Level: intermediate 2612 2613 Developer Notes: 2614 TSInterpolate() and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints. 2615 2616 .keywords: TS, set 2617 2618 .seealso: TSSetRetainStages(), TSSetPostStep() 2619 @*/ 2620 PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec U) 2621 { 2622 PetscErrorCode ierr; 2623 2624 PetscFunctionBegin; 2625 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2626 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 2627 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); 2628 if (!ts->ops->interpolate) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name); 2629 ierr = (*ts->ops->interpolate)(ts,t,U);CHKERRQ(ierr); 2630 PetscFunctionReturn(0); 2631 } 2632 2633 #undef __FUNCT__ 2634 #define __FUNCT__ "TSStep" 2635 /*@ 2636 TSStep - Steps one time step 2637 2638 Collective on TS 2639 2640 Input Parameter: 2641 . ts - the TS context obtained from TSCreate() 2642 2643 Level: intermediate 2644 2645 Notes: 2646 The hook set using TSSetPreStep() is called before each attempt to take the step. In general, the time step size may 2647 be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages. 2648 2649 This may over-step the final time provided in TSSetDuration() depending on the time-step used. TSSolve() interpolates to exactly the 2650 time provided in TSSetDuration(). One can use TSInterpolate() to determine an interpolated solution within the final timestep. 2651 2652 .keywords: TS, timestep, solve 2653 2654 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSInterpolate() 2655 @*/ 2656 PetscErrorCode TSStep(TS ts) 2657 { 2658 DM dm; 2659 PetscErrorCode ierr; 2660 static PetscBool cite = PETSC_FALSE; 2661 2662 PetscFunctionBegin; 2663 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2664 ierr = PetscCitationsRegister("@techreport{tspaper,\n" 2665 " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n" 2666 " author = {Shrirang Abhyankar and Jed Brown and Emil Constantinescu and Debojyoti Ghosh and Barry F. Smith},\n" 2667 " type = {Preprint},\n" 2668 " number = {ANL/MCS-P5061-0114},\n" 2669 " institution = {Argonne National Laboratory},\n" 2670 " year = {2014}\n}\n",&cite); 2671 2672 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 2673 ierr = TSSetUp(ts);CHKERRQ(ierr); 2674 2675 ts->reason = TS_CONVERGED_ITERATING; 2676 ts->ptime_prev = ts->ptime; 2677 ierr = DMSetOutputSequenceNumber(dm, ts->steps, ts->ptime);CHKERRQ(ierr); 2678 ierr = VecViewFromOptions(ts->vec_sol, ((PetscObject) ts)->prefix, "-ts_view_solution");CHKERRQ(ierr); 2679 2680 if (!ts->ops->step) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSStep not implemented for type '%s'",((PetscObject)ts)->type_name); 2681 ierr = PetscLogEventBegin(TS_Step,ts,0,0,0);CHKERRQ(ierr); 2682 if(ts->reverse_mode) { 2683 if(!ts->ops->stepadj) { 2684 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); 2685 }else { 2686 ierr = (*ts->ops->stepadj)(ts);CHKERRQ(ierr); 2687 } 2688 }else { 2689 ierr = (*ts->ops->step)(ts);CHKERRQ(ierr); 2690 } 2691 ierr = PetscLogEventEnd(TS_Step,ts,0,0,0);CHKERRQ(ierr); 2692 2693 ts->time_step_prev = ts->ptime - ts->ptime_prev; 2694 ierr = DMSetOutputSequenceNumber(dm, ts->steps, ts->ptime);CHKERRQ(ierr); 2695 2696 if (ts->reason < 0) { 2697 if (ts->errorifstepfailed) { 2698 if (ts->reason == TS_DIVERGED_NONLINEAR_SOLVE) { 2699 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]); 2700 } else if (ts->reason == TS_DIVERGED_STEP_REJECTED) { 2701 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]); 2702 } else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]); 2703 } 2704 } else if (!ts->reason) { 2705 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 2706 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 2707 } 2708 PetscFunctionReturn(0); 2709 } 2710 2711 #undef __FUNCT__ 2712 #define __FUNCT__ "TSEvaluateStep" 2713 /*@ 2714 TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy. 2715 2716 Collective on TS 2717 2718 Input Arguments: 2719 + ts - time stepping context 2720 . order - desired order of accuracy 2721 - done - whether the step was evaluated at this order (pass NULL to generate an error if not available) 2722 2723 Output Arguments: 2724 . U - state at the end of the current step 2725 2726 Level: advanced 2727 2728 Notes: 2729 This function cannot be called until all stages have been evaluated. 2730 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. 2731 2732 .seealso: TSStep(), TSAdapt 2733 @*/ 2734 PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec U,PetscBool *done) 2735 { 2736 PetscErrorCode ierr; 2737 2738 PetscFunctionBegin; 2739 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2740 PetscValidType(ts,1); 2741 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 2742 if (!ts->ops->evaluatestep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSEvaluateStep not implemented for type '%s'",((PetscObject)ts)->type_name); 2743 ierr = (*ts->ops->evaluatestep)(ts,order,U,done);CHKERRQ(ierr); 2744 PetscFunctionReturn(0); 2745 } 2746 2747 #undef __FUNCT__ 2748 #define __FUNCT__ "TSSolve" 2749 /*@ 2750 TSSolve - Steps the requested number of timesteps. 2751 2752 Collective on TS 2753 2754 Input Parameter: 2755 + ts - the TS context obtained from TSCreate() 2756 - u - the solution vector (can be null if TSSetSolution() was used, otherwise must contain the initial conditions) 2757 2758 Level: beginner 2759 2760 Notes: 2761 The final time returned by this function may be different from the time of the internally 2762 held state accessible by TSGetSolution() and TSGetTime() because the method may have 2763 stepped over the final time. 2764 2765 .keywords: TS, timestep, solve 2766 2767 .seealso: TSCreate(), TSSetSolution(), TSStep() 2768 @*/ 2769 PetscErrorCode TSSolve(TS ts,Vec u) 2770 { 2771 Vec solution; 2772 PetscErrorCode ierr; 2773 2774 PetscFunctionBegin; 2775 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2776 if (u) PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2777 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 */ 2778 PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2779 if (!ts->vec_sol || u == ts->vec_sol) { 2780 ierr = VecDuplicate(u,&solution);CHKERRQ(ierr); 2781 ierr = TSSetSolution(ts,solution);CHKERRQ(ierr); 2782 ierr = VecDestroy(&solution);CHKERRQ(ierr); /* grant ownership */ 2783 } 2784 ierr = VecCopy(u,ts->vec_sol);CHKERRQ(ierr); 2785 } else if (u) { 2786 ierr = TSSetSolution(ts,u);CHKERRQ(ierr); 2787 } 2788 ierr = TSSetUp(ts);CHKERRQ(ierr); /*compute adj coefficients if the reverse mode is on*/ 2789 /* reset time step and iteration counters */ 2790 ts->steps = 0; 2791 ts->ksp_its = 0; 2792 ts->snes_its = 0; 2793 ts->num_snes_failures = 0; 2794 ts->reject = 0; 2795 ts->reason = TS_CONVERGED_ITERATING; 2796 2797 ierr = TSViewFromOptions(ts,NULL,"-ts_view_pre");CHKERRQ(ierr); 2798 2799 if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */ 2800 ierr = (*ts->ops->solve)(ts);CHKERRQ(ierr); 2801 ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr); 2802 ts->solvetime = ts->ptime; 2803 } else { 2804 /* steps the requested number of timesteps. */ 2805 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 2806 else if (!ts->reverse_mode && ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 2807 while (!ts->reason) { 2808 if(!ts->reverse_mode) { 2809 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 2810 }else { 2811 ierr = TSMonitor(ts,ts->max_steps-ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 2812 } 2813 ierr = TSStep(ts);CHKERRQ(ierr); 2814 if (ts->event) { 2815 ierr = TSEventMonitor(ts);CHKERRQ(ierr); 2816 if (ts->event->status != TSEVENT_PROCESSING) { 2817 ierr = TSPostStep(ts);CHKERRQ(ierr); 2818 } 2819 } else { 2820 ierr = TSPostStep(ts);CHKERRQ(ierr); 2821 } 2822 } 2823 if (!ts->reverse_mode && ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) { 2824 ierr = TSInterpolate(ts,ts->max_time,u);CHKERRQ(ierr); 2825 ts->solvetime = ts->max_time; 2826 solution = u; 2827 } else { 2828 if (u) {ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr);} 2829 ts->solvetime = ts->ptime; 2830 solution = ts->vec_sol; 2831 } 2832 if(!ts->reverse_mode) { 2833 ierr = TSMonitor(ts,ts->steps,ts->solvetime,solution);CHKERRQ(ierr); 2834 } 2835 ierr = VecViewFromOptions(u, ((PetscObject) ts)->prefix, "-ts_view_solution");CHKERRQ(ierr); 2836 } 2837 2838 ierr = TSViewFromOptions(ts,NULL,"-ts_view");CHKERRQ(ierr); 2839 ierr = PetscObjectSAWsBlock((PetscObject)ts);CHKERRQ(ierr); 2840 PetscFunctionReturn(0); 2841 } 2842 2843 #undef __FUNCT__ 2844 #define __FUNCT__ "TSMonitor" 2845 /*@ 2846 TSMonitor - Runs all user-provided monitor routines set using TSMonitorSet() 2847 2848 Collective on TS 2849 2850 Input Parameters: 2851 + ts - time stepping context obtained from TSCreate() 2852 . step - step number that has just completed 2853 . ptime - model time of the state 2854 - u - state at the current model time 2855 2856 Notes: 2857 TSMonitor() is typically used within the time stepping implementations. 2858 Users might call this function when using the TSStep() interface instead of TSSolve(). 2859 2860 Level: advanced 2861 2862 .keywords: TS, timestep 2863 @*/ 2864 PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u) 2865 { 2866 PetscErrorCode ierr; 2867 PetscInt i,n = ts->numbermonitors; 2868 2869 PetscFunctionBegin; 2870 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2871 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 2872 ierr = VecLockPush(u);CHKERRQ(ierr); 2873 for (i=0; i<n; i++) { 2874 ierr = (*ts->monitor[i])(ts,step,ptime,u,ts->monitorcontext[i]);CHKERRQ(ierr); 2875 } 2876 ierr = VecLockPop(u);CHKERRQ(ierr); 2877 PetscFunctionReturn(0); 2878 } 2879 2880 /* ------------------------------------------------------------------------*/ 2881 #undef __FUNCT__ 2882 #define __FUNCT__ "TSMonitorLGCtxCreate" 2883 /*@C 2884 TSMonitorLGCtxCreate - Creates a line graph context for use with 2885 TS to monitor the solution process graphically in various ways 2886 2887 Collective on TS 2888 2889 Input Parameters: 2890 + host - the X display to open, or null for the local machine 2891 . label - the title to put in the title bar 2892 . x, y - the screen coordinates of the upper left coordinate of the window 2893 . m, n - the screen width and height in pixels 2894 - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time 2895 2896 Output Parameter: 2897 . ctx - the context 2898 2899 Options Database Key: 2900 + -ts_monitor_lg_timestep - automatically sets line graph monitor 2901 . -ts_monitor_lg_solution - 2902 . -ts_monitor_lg_error - 2903 . -ts_monitor_lg_ksp_iterations - 2904 . -ts_monitor_lg_snes_iterations - 2905 - -lg_indicate_data_points <true,false> - indicate the data points (at each time step) on the plot; default is true 2906 2907 Notes: 2908 Use TSMonitorLGCtxDestroy() to destroy. 2909 2910 Level: intermediate 2911 2912 .keywords: TS, monitor, line graph, residual, seealso 2913 2914 .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError() 2915 2916 @*/ 2917 PetscErrorCode TSMonitorLGCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorLGCtx *ctx) 2918 { 2919 PetscDraw win; 2920 PetscErrorCode ierr; 2921 2922 PetscFunctionBegin; 2923 ierr = PetscNew(ctx);CHKERRQ(ierr); 2924 ierr = PetscDrawCreate(comm,host,label,x,y,m,n,&win);CHKERRQ(ierr); 2925 ierr = PetscDrawSetFromOptions(win);CHKERRQ(ierr); 2926 ierr = PetscDrawLGCreate(win,1,&(*ctx)->lg);CHKERRQ(ierr); 2927 ierr = PetscLogObjectParent((PetscObject)(*ctx)->lg,(PetscObject)win);CHKERRQ(ierr); 2928 ierr = PetscDrawLGIndicateDataPoints((*ctx)->lg,PETSC_TRUE);CHKERRQ(ierr); 2929 ierr = PetscDrawLGSetFromOptions((*ctx)->lg);CHKERRQ(ierr); 2930 (*ctx)->howoften = howoften; 2931 PetscFunctionReturn(0); 2932 } 2933 2934 #undef __FUNCT__ 2935 #define __FUNCT__ "TSMonitorLGTimeStep" 2936 PetscErrorCode TSMonitorLGTimeStep(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx) 2937 { 2938 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 2939 PetscReal x = ptime,y; 2940 PetscErrorCode ierr; 2941 2942 PetscFunctionBegin; 2943 if (!step) { 2944 PetscDrawAxis axis; 2945 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 2946 ierr = PetscDrawAxisSetLabels(axis,"Timestep as function of time","Time","Time step");CHKERRQ(ierr); 2947 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 2948 ierr = PetscDrawLGIndicateDataPoints(ctx->lg,PETSC_TRUE);CHKERRQ(ierr); 2949 } 2950 ierr = TSGetTimeStep(ts,&y);CHKERRQ(ierr); 2951 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 2952 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 2953 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 2954 } 2955 PetscFunctionReturn(0); 2956 } 2957 2958 #undef __FUNCT__ 2959 #define __FUNCT__ "TSMonitorLGCtxDestroy" 2960 /*@C 2961 TSMonitorLGCtxDestroy - Destroys a line graph context that was created 2962 with TSMonitorLGCtxCreate(). 2963 2964 Collective on TSMonitorLGCtx 2965 2966 Input Parameter: 2967 . ctx - the monitor context 2968 2969 Level: intermediate 2970 2971 .keywords: TS, monitor, line graph, destroy 2972 2973 .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep(); 2974 @*/ 2975 PetscErrorCode TSMonitorLGCtxDestroy(TSMonitorLGCtx *ctx) 2976 { 2977 PetscDraw draw; 2978 PetscErrorCode ierr; 2979 2980 PetscFunctionBegin; 2981 ierr = PetscDrawLGGetDraw((*ctx)->lg,&draw);CHKERRQ(ierr); 2982 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 2983 ierr = PetscDrawLGDestroy(&(*ctx)->lg);CHKERRQ(ierr); 2984 ierr = PetscFree(*ctx);CHKERRQ(ierr); 2985 PetscFunctionReturn(0); 2986 } 2987 2988 #undef __FUNCT__ 2989 #define __FUNCT__ "TSGetTime" 2990 /*@ 2991 TSGetTime - Gets the time of the most recently completed step. 2992 2993 Not Collective 2994 2995 Input Parameter: 2996 . ts - the TS context obtained from TSCreate() 2997 2998 Output Parameter: 2999 . t - the current time 3000 3001 Level: beginner 3002 3003 Note: 3004 When called during time step evaluation (e.g. during residual evaluation or via hooks set using TSSetPreStep(), 3005 TSSetPreStage(), TSSetPostStage(), or TSSetPostStep()), the time is the time at the start of the step being evaluated. 3006 3007 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 3008 3009 .keywords: TS, get, time 3010 @*/ 3011 PetscErrorCode TSGetTime(TS ts,PetscReal *t) 3012 { 3013 PetscFunctionBegin; 3014 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3015 PetscValidRealPointer(t,2); 3016 *t = ts->ptime; 3017 PetscFunctionReturn(0); 3018 } 3019 3020 #undef __FUNCT__ 3021 #define __FUNCT__ "TSGetPrevTime" 3022 /*@ 3023 TSGetPrevTime - Gets the starting time of the previously completed step. 3024 3025 Not Collective 3026 3027 Input Parameter: 3028 . ts - the TS context obtained from TSCreate() 3029 3030 Output Parameter: 3031 . t - the previous time 3032 3033 Level: beginner 3034 3035 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 3036 3037 .keywords: TS, get, time 3038 @*/ 3039 PetscErrorCode TSGetPrevTime(TS ts,PetscReal *t) 3040 { 3041 PetscFunctionBegin; 3042 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3043 PetscValidRealPointer(t,2); 3044 *t = ts->ptime_prev; 3045 PetscFunctionReturn(0); 3046 } 3047 3048 #undef __FUNCT__ 3049 #define __FUNCT__ "TSSetTime" 3050 /*@ 3051 TSSetTime - Allows one to reset the time. 3052 3053 Logically Collective on TS 3054 3055 Input Parameters: 3056 + ts - the TS context obtained from TSCreate() 3057 - time - the time 3058 3059 Level: intermediate 3060 3061 .seealso: TSGetTime(), TSSetDuration() 3062 3063 .keywords: TS, set, time 3064 @*/ 3065 PetscErrorCode TSSetTime(TS ts, PetscReal t) 3066 { 3067 PetscFunctionBegin; 3068 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3069 PetscValidLogicalCollectiveReal(ts,t,2); 3070 ts->ptime = t; 3071 PetscFunctionReturn(0); 3072 } 3073 3074 #undef __FUNCT__ 3075 #define __FUNCT__ "TSSetOptionsPrefix" 3076 /*@C 3077 TSSetOptionsPrefix - Sets the prefix used for searching for all 3078 TS options in the database. 3079 3080 Logically Collective on TS 3081 3082 Input Parameter: 3083 + ts - The TS context 3084 - prefix - The prefix to prepend to all option names 3085 3086 Notes: 3087 A hyphen (-) must NOT be given at the beginning of the prefix name. 3088 The first character of all runtime options is AUTOMATICALLY the 3089 hyphen. 3090 3091 Level: advanced 3092 3093 .keywords: TS, set, options, prefix, database 3094 3095 .seealso: TSSetFromOptions() 3096 3097 @*/ 3098 PetscErrorCode TSSetOptionsPrefix(TS ts,const char prefix[]) 3099 { 3100 PetscErrorCode ierr; 3101 SNES snes; 3102 3103 PetscFunctionBegin; 3104 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3105 ierr = PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 3106 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3107 ierr = SNESSetOptionsPrefix(snes,prefix);CHKERRQ(ierr); 3108 PetscFunctionReturn(0); 3109 } 3110 3111 3112 #undef __FUNCT__ 3113 #define __FUNCT__ "TSAppendOptionsPrefix" 3114 /*@C 3115 TSAppendOptionsPrefix - Appends to the prefix used for searching for all 3116 TS options in the database. 3117 3118 Logically Collective on TS 3119 3120 Input Parameter: 3121 + ts - The TS context 3122 - prefix - The prefix to prepend to all option names 3123 3124 Notes: 3125 A hyphen (-) must NOT be given at the beginning of the prefix name. 3126 The first character of all runtime options is AUTOMATICALLY the 3127 hyphen. 3128 3129 Level: advanced 3130 3131 .keywords: TS, append, options, prefix, database 3132 3133 .seealso: TSGetOptionsPrefix() 3134 3135 @*/ 3136 PetscErrorCode TSAppendOptionsPrefix(TS ts,const char prefix[]) 3137 { 3138 PetscErrorCode ierr; 3139 SNES snes; 3140 3141 PetscFunctionBegin; 3142 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3143 ierr = PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 3144 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3145 ierr = SNESAppendOptionsPrefix(snes,prefix);CHKERRQ(ierr); 3146 PetscFunctionReturn(0); 3147 } 3148 3149 #undef __FUNCT__ 3150 #define __FUNCT__ "TSGetOptionsPrefix" 3151 /*@C 3152 TSGetOptionsPrefix - Sets the prefix used for searching for all 3153 TS options in the database. 3154 3155 Not Collective 3156 3157 Input Parameter: 3158 . ts - The TS context 3159 3160 Output Parameter: 3161 . prefix - A pointer to the prefix string used 3162 3163 Notes: On the fortran side, the user should pass in a string 'prifix' of 3164 sufficient length to hold the prefix. 3165 3166 Level: intermediate 3167 3168 .keywords: TS, get, options, prefix, database 3169 3170 .seealso: TSAppendOptionsPrefix() 3171 @*/ 3172 PetscErrorCode TSGetOptionsPrefix(TS ts,const char *prefix[]) 3173 { 3174 PetscErrorCode ierr; 3175 3176 PetscFunctionBegin; 3177 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3178 PetscValidPointer(prefix,2); 3179 ierr = PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 3180 PetscFunctionReturn(0); 3181 } 3182 3183 #undef __FUNCT__ 3184 #define __FUNCT__ "TSGetRHSJacobian" 3185 /*@C 3186 TSGetRHSJacobian - Returns the Jacobian J at the present timestep. 3187 3188 Not Collective, but parallel objects are returned if TS is parallel 3189 3190 Input Parameter: 3191 . ts - The TS context obtained from TSCreate() 3192 3193 Output Parameters: 3194 + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or NULL) 3195 . Pmat - The matrix from which the preconditioner is constructed, usually the same as Amat (or NULL) 3196 . func - Function to compute the Jacobian of the RHS (or NULL) 3197 - ctx - User-defined context for Jacobian evaluation routine (or NULL) 3198 3199 Notes: You can pass in NULL for any return argument you do not need. 3200 3201 Level: intermediate 3202 3203 .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 3204 3205 .keywords: TS, timestep, get, matrix, Jacobian 3206 @*/ 3207 PetscErrorCode TSGetRHSJacobian(TS ts,Mat *Amat,Mat *Pmat,TSRHSJacobian *func,void **ctx) 3208 { 3209 PetscErrorCode ierr; 3210 SNES snes; 3211 DM dm; 3212 3213 PetscFunctionBegin; 3214 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3215 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 3216 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 3217 ierr = DMTSGetRHSJacobian(dm,func,ctx);CHKERRQ(ierr); 3218 PetscFunctionReturn(0); 3219 } 3220 3221 #undef __FUNCT__ 3222 #define __FUNCT__ "TSGetIJacobian" 3223 /*@C 3224 TSGetIJacobian - Returns the implicit Jacobian at the present timestep. 3225 3226 Not Collective, but parallel objects are returned if TS is parallel 3227 3228 Input Parameter: 3229 . ts - The TS context obtained from TSCreate() 3230 3231 Output Parameters: 3232 + Amat - The (approximate) Jacobian of F(t,U,U_t) 3233 . Pmat - The matrix from which the preconditioner is constructed, often the same as Amat 3234 . f - The function to compute the matrices 3235 - ctx - User-defined context for Jacobian evaluation routine 3236 3237 Notes: You can pass in NULL for any return argument you do not need. 3238 3239 Level: advanced 3240 3241 .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 3242 3243 .keywords: TS, timestep, get, matrix, Jacobian 3244 @*/ 3245 PetscErrorCode TSGetIJacobian(TS ts,Mat *Amat,Mat *Pmat,TSIJacobian *f,void **ctx) 3246 { 3247 PetscErrorCode ierr; 3248 SNES snes; 3249 DM dm; 3250 3251 PetscFunctionBegin; 3252 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3253 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 3254 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 3255 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 3256 ierr = DMTSGetIJacobian(dm,f,ctx);CHKERRQ(ierr); 3257 PetscFunctionReturn(0); 3258 } 3259 3260 3261 #undef __FUNCT__ 3262 #define __FUNCT__ "TSMonitorDrawSolution" 3263 /*@C 3264 TSMonitorDrawSolution - Monitors progress of the TS solvers by calling 3265 VecView() for the solution at each timestep 3266 3267 Collective on TS 3268 3269 Input Parameters: 3270 + ts - the TS context 3271 . step - current time-step 3272 . ptime - current time 3273 - dummy - either a viewer or NULL 3274 3275 Options Database: 3276 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 3277 3278 Notes: the initial solution and current solution are not displayed with a common axis scaling so generally the option -ts_monitor_draw_solution_initial 3279 will look bad 3280 3281 Level: intermediate 3282 3283 .keywords: TS, vector, monitor, view 3284 3285 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3286 @*/ 3287 PetscErrorCode TSMonitorDrawSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3288 { 3289 PetscErrorCode ierr; 3290 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy; 3291 PetscDraw draw; 3292 3293 PetscFunctionBegin; 3294 if (!step && ictx->showinitial) { 3295 if (!ictx->initialsolution) { 3296 ierr = VecDuplicate(u,&ictx->initialsolution);CHKERRQ(ierr); 3297 } 3298 ierr = VecCopy(u,ictx->initialsolution);CHKERRQ(ierr); 3299 } 3300 if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) PetscFunctionReturn(0); 3301 3302 if (ictx->showinitial) { 3303 PetscReal pause; 3304 ierr = PetscViewerDrawGetPause(ictx->viewer,&pause);CHKERRQ(ierr); 3305 ierr = PetscViewerDrawSetPause(ictx->viewer,0.0);CHKERRQ(ierr); 3306 ierr = VecView(ictx->initialsolution,ictx->viewer);CHKERRQ(ierr); 3307 ierr = PetscViewerDrawSetPause(ictx->viewer,pause);CHKERRQ(ierr); 3308 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);CHKERRQ(ierr); 3309 } 3310 ierr = VecView(u,ictx->viewer);CHKERRQ(ierr); 3311 if (ictx->showtimestepandtime) { 3312 PetscReal xl,yl,xr,yr,tw,w,h; 3313 char time[32]; 3314 size_t len; 3315 3316 ierr = PetscViewerDrawGetDraw(ictx->viewer,0,&draw);CHKERRQ(ierr); 3317 ierr = PetscSNPrintf(time,32,"Timestep %d Time %f",(int)step,(double)ptime);CHKERRQ(ierr); 3318 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 3319 ierr = PetscStrlen(time,&len);CHKERRQ(ierr); 3320 ierr = PetscDrawStringGetSize(draw,&tw,NULL);CHKERRQ(ierr); 3321 w = xl + .5*(xr - xl) - .5*len*tw; 3322 h = yl + .95*(yr - yl); 3323 ierr = PetscDrawString(draw,w,h,PETSC_DRAW_BLACK,time);CHKERRQ(ierr); 3324 ierr = PetscDrawFlush(draw);CHKERRQ(ierr); 3325 } 3326 3327 if (ictx->showinitial) { 3328 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);CHKERRQ(ierr); 3329 } 3330 PetscFunctionReturn(0); 3331 } 3332 3333 #undef __FUNCT__ 3334 #define __FUNCT__ "TSMonitorDrawSolutionPhase" 3335 /*@C 3336 TSMonitorDrawSolutionPhase - Monitors progress of the TS solvers by plotting the solution as a phase diagram 3337 3338 Collective on TS 3339 3340 Input Parameters: 3341 + ts - the TS context 3342 . step - current time-step 3343 . ptime - current time 3344 - dummy - either a viewer or NULL 3345 3346 Level: intermediate 3347 3348 .keywords: TS, vector, monitor, view 3349 3350 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3351 @*/ 3352 PetscErrorCode TSMonitorDrawSolutionPhase(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3353 { 3354 PetscErrorCode ierr; 3355 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy; 3356 PetscDraw draw; 3357 MPI_Comm comm; 3358 PetscInt n; 3359 PetscMPIInt size; 3360 PetscReal xl,yl,xr,yr,tw,w,h; 3361 char time[32]; 3362 size_t len; 3363 const PetscScalar *U; 3364 3365 PetscFunctionBegin; 3366 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 3367 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3368 if (size != 1) SETERRQ(comm,PETSC_ERR_SUP,"Only allowed for sequential runs"); 3369 ierr = VecGetSize(u,&n);CHKERRQ(ierr); 3370 if (n != 2) SETERRQ(comm,PETSC_ERR_SUP,"Only for ODEs with two unknowns"); 3371 3372 ierr = PetscViewerDrawGetDraw(ictx->viewer,0,&draw);CHKERRQ(ierr); 3373 3374 ierr = VecGetArrayRead(u,&U);CHKERRQ(ierr); 3375 ierr = PetscDrawAxisGetLimits(ictx->axis,&xl,&xr,&yl,&yr);CHKERRQ(ierr); 3376 if ((PetscRealPart(U[0]) < xl) || (PetscRealPart(U[1]) < yl) || (PetscRealPart(U[0]) > xr) || (PetscRealPart(U[1]) > yr)) { 3377 ierr = VecRestoreArrayRead(u,&U);CHKERRQ(ierr); 3378 PetscFunctionReturn(0); 3379 } 3380 if (!step) ictx->color++; 3381 ierr = PetscDrawPoint(draw,PetscRealPart(U[0]),PetscRealPart(U[1]),ictx->color);CHKERRQ(ierr); 3382 ierr = VecRestoreArrayRead(u,&U);CHKERRQ(ierr); 3383 3384 if (ictx->showtimestepandtime) { 3385 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 3386 ierr = PetscSNPrintf(time,32,"Timestep %d Time %f",(int)step,(double)ptime);CHKERRQ(ierr); 3387 ierr = PetscStrlen(time,&len);CHKERRQ(ierr); 3388 ierr = PetscDrawStringGetSize(draw,&tw,NULL);CHKERRQ(ierr); 3389 w = xl + .5*(xr - xl) - .5*len*tw; 3390 h = yl + .95*(yr - yl); 3391 ierr = PetscDrawString(draw,w,h,PETSC_DRAW_BLACK,time);CHKERRQ(ierr); 3392 } 3393 ierr = PetscDrawFlush(draw);CHKERRQ(ierr); 3394 PetscFunctionReturn(0); 3395 } 3396 3397 3398 #undef __FUNCT__ 3399 #define __FUNCT__ "TSMonitorDrawCtxDestroy" 3400 /*@C 3401 TSMonitorDrawCtxDestroy - Destroys the monitor context for TSMonitorDrawSolution() 3402 3403 Collective on TS 3404 3405 Input Parameters: 3406 . ctx - the monitor context 3407 3408 Level: intermediate 3409 3410 .keywords: TS, vector, monitor, view 3411 3412 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawSolution(), TSMonitorDrawError() 3413 @*/ 3414 PetscErrorCode TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx) 3415 { 3416 PetscErrorCode ierr; 3417 3418 PetscFunctionBegin; 3419 ierr = PetscDrawAxisDestroy(&(*ictx)->axis);CHKERRQ(ierr); 3420 ierr = PetscViewerDestroy(&(*ictx)->viewer);CHKERRQ(ierr); 3421 ierr = VecDestroy(&(*ictx)->initialsolution);CHKERRQ(ierr); 3422 ierr = PetscFree(*ictx);CHKERRQ(ierr); 3423 PetscFunctionReturn(0); 3424 } 3425 3426 #undef __FUNCT__ 3427 #define __FUNCT__ "TSMonitorDrawCtxCreate" 3428 /*@C 3429 TSMonitorDrawCtxCreate - Creates the monitor context for TSMonitorDrawCtx 3430 3431 Collective on TS 3432 3433 Input Parameter: 3434 . ts - time-step context 3435 3436 Output Patameter: 3437 . ctx - the monitor context 3438 3439 Options Database: 3440 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 3441 3442 Level: intermediate 3443 3444 .keywords: TS, vector, monitor, view 3445 3446 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawCtx() 3447 @*/ 3448 PetscErrorCode TSMonitorDrawCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorDrawCtx *ctx) 3449 { 3450 PetscErrorCode ierr; 3451 3452 PetscFunctionBegin; 3453 ierr = PetscNew(ctx);CHKERRQ(ierr); 3454 ierr = PetscViewerDrawOpen(comm,host,label,x,y,m,n,&(*ctx)->viewer);CHKERRQ(ierr); 3455 ierr = PetscViewerSetFromOptions((*ctx)->viewer);CHKERRQ(ierr); 3456 3457 (*ctx)->howoften = howoften; 3458 (*ctx)->showinitial = PETSC_FALSE; 3459 ierr = PetscOptionsGetBool(NULL,"-ts_monitor_draw_solution_initial",&(*ctx)->showinitial,NULL);CHKERRQ(ierr); 3460 3461 (*ctx)->showtimestepandtime = PETSC_FALSE; 3462 ierr = PetscOptionsGetBool(NULL,"-ts_monitor_draw_solution_show_time",&(*ctx)->showtimestepandtime,NULL);CHKERRQ(ierr); 3463 (*ctx)->color = PETSC_DRAW_WHITE; 3464 PetscFunctionReturn(0); 3465 } 3466 3467 #undef __FUNCT__ 3468 #define __FUNCT__ "TSMonitorDrawError" 3469 /*@C 3470 TSMonitorDrawError - Monitors progress of the TS solvers by calling 3471 VecView() for the error at each timestep 3472 3473 Collective on TS 3474 3475 Input Parameters: 3476 + ts - the TS context 3477 . step - current time-step 3478 . ptime - current time 3479 - dummy - either a viewer or NULL 3480 3481 Level: intermediate 3482 3483 .keywords: TS, vector, monitor, view 3484 3485 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3486 @*/ 3487 PetscErrorCode TSMonitorDrawError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3488 { 3489 PetscErrorCode ierr; 3490 TSMonitorDrawCtx ctx = (TSMonitorDrawCtx)dummy; 3491 PetscViewer viewer = ctx->viewer; 3492 Vec work; 3493 3494 PetscFunctionBegin; 3495 if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) PetscFunctionReturn(0); 3496 ierr = VecDuplicate(u,&work);CHKERRQ(ierr); 3497 ierr = TSComputeSolutionFunction(ts,ptime,work);CHKERRQ(ierr); 3498 ierr = VecAXPY(work,-1.0,u);CHKERRQ(ierr); 3499 ierr = VecView(work,viewer);CHKERRQ(ierr); 3500 ierr = VecDestroy(&work);CHKERRQ(ierr); 3501 PetscFunctionReturn(0); 3502 } 3503 3504 #include <petsc-private/dmimpl.h> 3505 #undef __FUNCT__ 3506 #define __FUNCT__ "TSSetDM" 3507 /*@ 3508 TSSetDM - Sets the DM that may be used by some preconditioners 3509 3510 Logically Collective on TS and DM 3511 3512 Input Parameters: 3513 + ts - the preconditioner context 3514 - dm - the dm 3515 3516 Level: intermediate 3517 3518 3519 .seealso: TSGetDM(), SNESSetDM(), SNESGetDM() 3520 @*/ 3521 PetscErrorCode TSSetDM(TS ts,DM dm) 3522 { 3523 PetscErrorCode ierr; 3524 SNES snes; 3525 DMTS tsdm; 3526 3527 PetscFunctionBegin; 3528 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3529 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 3530 if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */ 3531 if (ts->dm->dmts && !dm->dmts) { 3532 ierr = DMCopyDMTS(ts->dm,dm);CHKERRQ(ierr); 3533 ierr = DMGetDMTS(ts->dm,&tsdm);CHKERRQ(ierr); 3534 if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */ 3535 tsdm->originaldm = dm; 3536 } 3537 } 3538 ierr = DMDestroy(&ts->dm);CHKERRQ(ierr); 3539 } 3540 ts->dm = dm; 3541 3542 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3543 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 3544 PetscFunctionReturn(0); 3545 } 3546 3547 #undef __FUNCT__ 3548 #define __FUNCT__ "TSGetDM" 3549 /*@ 3550 TSGetDM - Gets the DM that may be used by some preconditioners 3551 3552 Not Collective 3553 3554 Input Parameter: 3555 . ts - the preconditioner context 3556 3557 Output Parameter: 3558 . dm - the dm 3559 3560 Level: intermediate 3561 3562 3563 .seealso: TSSetDM(), SNESSetDM(), SNESGetDM() 3564 @*/ 3565 PetscErrorCode TSGetDM(TS ts,DM *dm) 3566 { 3567 PetscErrorCode ierr; 3568 3569 PetscFunctionBegin; 3570 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3571 if (!ts->dm) { 3572 ierr = DMShellCreate(PetscObjectComm((PetscObject)ts),&ts->dm);CHKERRQ(ierr); 3573 if (ts->snes) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 3574 } 3575 *dm = ts->dm; 3576 PetscFunctionReturn(0); 3577 } 3578 3579 #undef __FUNCT__ 3580 #define __FUNCT__ "SNESTSFormFunction" 3581 /*@ 3582 SNESTSFormFunction - Function to evaluate nonlinear residual 3583 3584 Logically Collective on SNES 3585 3586 Input Parameter: 3587 + snes - nonlinear solver 3588 . U - the current state at which to evaluate the residual 3589 - ctx - user context, must be a TS 3590 3591 Output Parameter: 3592 . F - the nonlinear residual 3593 3594 Notes: 3595 This function is not normally called by users and is automatically registered with the SNES used by TS. 3596 It is most frequently passed to MatFDColoringSetFunction(). 3597 3598 Level: advanced 3599 3600 .seealso: SNESSetFunction(), MatFDColoringSetFunction() 3601 @*/ 3602 PetscErrorCode SNESTSFormFunction(SNES snes,Vec U,Vec F,void *ctx) 3603 { 3604 TS ts = (TS)ctx; 3605 PetscErrorCode ierr; 3606 3607 PetscFunctionBegin; 3608 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3609 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 3610 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 3611 PetscValidHeaderSpecific(ts,TS_CLASSID,4); 3612 ierr = (ts->ops->snesfunction)(snes,U,F,ts);CHKERRQ(ierr); 3613 PetscFunctionReturn(0); 3614 } 3615 3616 #undef __FUNCT__ 3617 #define __FUNCT__ "SNESTSFormJacobian" 3618 /*@ 3619 SNESTSFormJacobian - Function to evaluate the Jacobian 3620 3621 Collective on SNES 3622 3623 Input Parameter: 3624 + snes - nonlinear solver 3625 . U - the current state at which to evaluate the residual 3626 - ctx - user context, must be a TS 3627 3628 Output Parameter: 3629 + A - the Jacobian 3630 . B - the preconditioning matrix (may be the same as A) 3631 - flag - indicates any structure change in the matrix 3632 3633 Notes: 3634 This function is not normally called by users and is automatically registered with the SNES used by TS. 3635 3636 Level: developer 3637 3638 .seealso: SNESSetJacobian() 3639 @*/ 3640 PetscErrorCode SNESTSFormJacobian(SNES snes,Vec U,Mat A,Mat B,void *ctx) 3641 { 3642 TS ts = (TS)ctx; 3643 PetscErrorCode ierr; 3644 3645 PetscFunctionBegin; 3646 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3647 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 3648 PetscValidPointer(A,3); 3649 PetscValidHeaderSpecific(A,MAT_CLASSID,3); 3650 PetscValidPointer(B,4); 3651 PetscValidHeaderSpecific(B,MAT_CLASSID,4); 3652 PetscValidHeaderSpecific(ts,TS_CLASSID,6); 3653 ierr = (ts->ops->snesjacobian)(snes,U,A,B,ts);CHKERRQ(ierr); 3654 PetscFunctionReturn(0); 3655 } 3656 3657 #undef __FUNCT__ 3658 #define __FUNCT__ "TSComputeRHSFunctionLinear" 3659 /*@C 3660 TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems only 3661 3662 Collective on TS 3663 3664 Input Arguments: 3665 + ts - time stepping context 3666 . t - time at which to evaluate 3667 . U - state at which to evaluate 3668 - ctx - context 3669 3670 Output Arguments: 3671 . F - right hand side 3672 3673 Level: intermediate 3674 3675 Notes: 3676 This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems. 3677 The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian(). 3678 3679 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant() 3680 @*/ 3681 PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec U,Vec F,void *ctx) 3682 { 3683 PetscErrorCode ierr; 3684 Mat Arhs,Brhs; 3685 3686 PetscFunctionBegin; 3687 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 3688 ierr = TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);CHKERRQ(ierr); 3689 ierr = MatMult(Arhs,U,F);CHKERRQ(ierr); 3690 PetscFunctionReturn(0); 3691 } 3692 3693 #undef __FUNCT__ 3694 #define __FUNCT__ "TSComputeRHSJacobianConstant" 3695 /*@C 3696 TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent. 3697 3698 Collective on TS 3699 3700 Input Arguments: 3701 + ts - time stepping context 3702 . t - time at which to evaluate 3703 . U - state at which to evaluate 3704 - ctx - context 3705 3706 Output Arguments: 3707 + A - pointer to operator 3708 . B - pointer to preconditioning matrix 3709 - flg - matrix structure flag 3710 3711 Level: intermediate 3712 3713 Notes: 3714 This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems. 3715 3716 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear() 3717 @*/ 3718 PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx) 3719 { 3720 PetscFunctionBegin; 3721 PetscFunctionReturn(0); 3722 } 3723 3724 #undef __FUNCT__ 3725 #define __FUNCT__ "TSComputeIFunctionLinear" 3726 /*@C 3727 TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only 3728 3729 Collective on TS 3730 3731 Input Arguments: 3732 + ts - time stepping context 3733 . t - time at which to evaluate 3734 . U - state at which to evaluate 3735 . Udot - time derivative of state vector 3736 - ctx - context 3737 3738 Output Arguments: 3739 . F - left hand side 3740 3741 Level: intermediate 3742 3743 Notes: 3744 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 3745 user is required to write their own TSComputeIFunction. 3746 This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems. 3747 The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian(). 3748 3749 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant() 3750 @*/ 3751 PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx) 3752 { 3753 PetscErrorCode ierr; 3754 Mat A,B; 3755 3756 PetscFunctionBegin; 3757 ierr = TSGetIJacobian(ts,&A,&B,NULL,NULL);CHKERRQ(ierr); 3758 ierr = TSComputeIJacobian(ts,t,U,Udot,1.0,A,B,PETSC_TRUE);CHKERRQ(ierr); 3759 ierr = MatMult(A,Udot,F);CHKERRQ(ierr); 3760 PetscFunctionReturn(0); 3761 } 3762 3763 #undef __FUNCT__ 3764 #define __FUNCT__ "TSComputeIJacobianConstant" 3765 /*@C 3766 TSComputeIJacobianConstant - Reuses a time-independent for a semi-implicit DAE or ODE 3767 3768 Collective on TS 3769 3770 Input Arguments: 3771 + ts - time stepping context 3772 . t - time at which to evaluate 3773 . U - state at which to evaluate 3774 . Udot - time derivative of state vector 3775 . shift - shift to apply 3776 - ctx - context 3777 3778 Output Arguments: 3779 + A - pointer to operator 3780 . B - pointer to preconditioning matrix 3781 - flg - matrix structure flag 3782 3783 Level: advanced 3784 3785 Notes: 3786 This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems. 3787 3788 It is only appropriate for problems of the form 3789 3790 $ M Udot = F(U,t) 3791 3792 where M is constant and F is non-stiff. The user must pass M to TSSetIJacobian(). The current implementation only 3793 works with IMEX time integration methods such as TSROSW and TSARKIMEX, since there is no support for de-constructing 3794 an implicit operator of the form 3795 3796 $ shift*M + J 3797 3798 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 3799 a copy of M or reassemble it when requested. 3800 3801 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear() 3802 @*/ 3803 PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat A,Mat B,void *ctx) 3804 { 3805 PetscErrorCode ierr; 3806 3807 PetscFunctionBegin; 3808 ierr = MatScale(A, shift / ts->ijacobian.shift);CHKERRQ(ierr); 3809 ts->ijacobian.shift = shift; 3810 PetscFunctionReturn(0); 3811 } 3812 3813 #undef __FUNCT__ 3814 #define __FUNCT__ "TSGetEquationType" 3815 /*@ 3816 TSGetEquationType - Gets the type of the equation that TS is solving. 3817 3818 Not Collective 3819 3820 Input Parameter: 3821 . ts - the TS context 3822 3823 Output Parameter: 3824 . equation_type - see TSEquationType 3825 3826 Level: beginner 3827 3828 .keywords: TS, equation type 3829 3830 .seealso: TSSetEquationType(), TSEquationType 3831 @*/ 3832 PetscErrorCode TSGetEquationType(TS ts,TSEquationType *equation_type) 3833 { 3834 PetscFunctionBegin; 3835 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3836 PetscValidPointer(equation_type,2); 3837 *equation_type = ts->equation_type; 3838 PetscFunctionReturn(0); 3839 } 3840 3841 #undef __FUNCT__ 3842 #define __FUNCT__ "TSSetEquationType" 3843 /*@ 3844 TSSetEquationType - Sets the type of the equation that TS is solving. 3845 3846 Not Collective 3847 3848 Input Parameter: 3849 + ts - the TS context 3850 . equation_type - see TSEquationType 3851 3852 Level: advanced 3853 3854 .keywords: TS, equation type 3855 3856 .seealso: TSGetEquationType(), TSEquationType 3857 @*/ 3858 PetscErrorCode TSSetEquationType(TS ts,TSEquationType equation_type) 3859 { 3860 PetscFunctionBegin; 3861 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3862 ts->equation_type = equation_type; 3863 PetscFunctionReturn(0); 3864 } 3865 3866 #undef __FUNCT__ 3867 #define __FUNCT__ "TSGetConvergedReason" 3868 /*@ 3869 TSGetConvergedReason - Gets the reason the TS iteration was stopped. 3870 3871 Not Collective 3872 3873 Input Parameter: 3874 . ts - the TS context 3875 3876 Output Parameter: 3877 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 3878 manual pages for the individual convergence tests for complete lists 3879 3880 Level: beginner 3881 3882 Notes: 3883 Can only be called after the call to TSSolve() is complete. 3884 3885 .keywords: TS, nonlinear, set, convergence, test 3886 3887 .seealso: TSSetConvergenceTest(), TSConvergedReason 3888 @*/ 3889 PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason) 3890 { 3891 PetscFunctionBegin; 3892 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3893 PetscValidPointer(reason,2); 3894 *reason = ts->reason; 3895 PetscFunctionReturn(0); 3896 } 3897 3898 #undef __FUNCT__ 3899 #define __FUNCT__ "TSSetConvergedReason" 3900 /*@ 3901 TSSetConvergedReason - Sets the reason for handling the convergence of TSSolve. 3902 3903 Not Collective 3904 3905 Input Parameter: 3906 + ts - the TS context 3907 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 3908 manual pages for the individual convergence tests for complete lists 3909 3910 Level: advanced 3911 3912 Notes: 3913 Can only be called during TSSolve() is active. 3914 3915 .keywords: TS, nonlinear, set, convergence, test 3916 3917 .seealso: TSConvergedReason 3918 @*/ 3919 PetscErrorCode TSSetConvergedReason(TS ts,TSConvergedReason reason) 3920 { 3921 PetscFunctionBegin; 3922 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3923 ts->reason = reason; 3924 PetscFunctionReturn(0); 3925 } 3926 3927 #undef __FUNCT__ 3928 #define __FUNCT__ "TSGetSolveTime" 3929 /*@ 3930 TSGetSolveTime - Gets the time after a call to TSSolve() 3931 3932 Not Collective 3933 3934 Input Parameter: 3935 . ts - the TS context 3936 3937 Output Parameter: 3938 . ftime - the final time. This time should correspond to the final time set with TSSetDuration() 3939 3940 Level: beginner 3941 3942 Notes: 3943 Can only be called after the call to TSSolve() is complete. 3944 3945 .keywords: TS, nonlinear, set, convergence, test 3946 3947 .seealso: TSSetConvergenceTest(), TSConvergedReason 3948 @*/ 3949 PetscErrorCode TSGetSolveTime(TS ts,PetscReal *ftime) 3950 { 3951 PetscFunctionBegin; 3952 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3953 PetscValidPointer(ftime,2); 3954 *ftime = ts->solvetime; 3955 PetscFunctionReturn(0); 3956 } 3957 3958 #undef __FUNCT__ 3959 #define __FUNCT__ "TSGetSNESIterations" 3960 /*@ 3961 TSGetSNESIterations - Gets the total number of nonlinear iterations 3962 used by the time integrator. 3963 3964 Not Collective 3965 3966 Input Parameter: 3967 . ts - TS context 3968 3969 Output Parameter: 3970 . nits - number of nonlinear iterations 3971 3972 Notes: 3973 This counter is reset to zero for each successive call to TSSolve(). 3974 3975 Level: intermediate 3976 3977 .keywords: TS, get, number, nonlinear, iterations 3978 3979 .seealso: TSGetKSPIterations() 3980 @*/ 3981 PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits) 3982 { 3983 PetscFunctionBegin; 3984 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3985 PetscValidIntPointer(nits,2); 3986 *nits = ts->snes_its; 3987 PetscFunctionReturn(0); 3988 } 3989 3990 #undef __FUNCT__ 3991 #define __FUNCT__ "TSGetKSPIterations" 3992 /*@ 3993 TSGetKSPIterations - Gets the total number of linear iterations 3994 used by the time integrator. 3995 3996 Not Collective 3997 3998 Input Parameter: 3999 . ts - TS context 4000 4001 Output Parameter: 4002 . lits - number of linear iterations 4003 4004 Notes: 4005 This counter is reset to zero for each successive call to TSSolve(). 4006 4007 Level: intermediate 4008 4009 .keywords: TS, get, number, linear, iterations 4010 4011 .seealso: TSGetSNESIterations(), SNESGetKSPIterations() 4012 @*/ 4013 PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits) 4014 { 4015 PetscFunctionBegin; 4016 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4017 PetscValidIntPointer(lits,2); 4018 *lits = ts->ksp_its; 4019 PetscFunctionReturn(0); 4020 } 4021 4022 #undef __FUNCT__ 4023 #define __FUNCT__ "TSGetStepRejections" 4024 /*@ 4025 TSGetStepRejections - Gets the total number of rejected steps. 4026 4027 Not Collective 4028 4029 Input Parameter: 4030 . ts - TS context 4031 4032 Output Parameter: 4033 . rejects - number of steps rejected 4034 4035 Notes: 4036 This counter is reset to zero for each successive call to TSSolve(). 4037 4038 Level: intermediate 4039 4040 .keywords: TS, get, number 4041 4042 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails() 4043 @*/ 4044 PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects) 4045 { 4046 PetscFunctionBegin; 4047 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4048 PetscValidIntPointer(rejects,2); 4049 *rejects = ts->reject; 4050 PetscFunctionReturn(0); 4051 } 4052 4053 #undef __FUNCT__ 4054 #define __FUNCT__ "TSGetSNESFailures" 4055 /*@ 4056 TSGetSNESFailures - Gets the total number of failed SNES solves 4057 4058 Not Collective 4059 4060 Input Parameter: 4061 . ts - TS context 4062 4063 Output Parameter: 4064 . fails - number of failed nonlinear solves 4065 4066 Notes: 4067 This counter is reset to zero for each successive call to TSSolve(). 4068 4069 Level: intermediate 4070 4071 .keywords: TS, get, number 4072 4073 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures() 4074 @*/ 4075 PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails) 4076 { 4077 PetscFunctionBegin; 4078 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4079 PetscValidIntPointer(fails,2); 4080 *fails = ts->num_snes_failures; 4081 PetscFunctionReturn(0); 4082 } 4083 4084 #undef __FUNCT__ 4085 #define __FUNCT__ "TSSetMaxStepRejections" 4086 /*@ 4087 TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails 4088 4089 Not Collective 4090 4091 Input Parameter: 4092 + ts - TS context 4093 - rejects - maximum number of rejected steps, pass -1 for unlimited 4094 4095 Notes: 4096 The counter is reset to zero for each step 4097 4098 Options Database Key: 4099 . -ts_max_reject - Maximum number of step rejections before a step fails 4100 4101 Level: intermediate 4102 4103 .keywords: TS, set, maximum, number 4104 4105 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 4106 @*/ 4107 PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects) 4108 { 4109 PetscFunctionBegin; 4110 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4111 ts->max_reject = rejects; 4112 PetscFunctionReturn(0); 4113 } 4114 4115 #undef __FUNCT__ 4116 #define __FUNCT__ "TSSetMaxSNESFailures" 4117 /*@ 4118 TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves 4119 4120 Not Collective 4121 4122 Input Parameter: 4123 + ts - TS context 4124 - fails - maximum number of failed nonlinear solves, pass -1 for unlimited 4125 4126 Notes: 4127 The counter is reset to zero for each successive call to TSSolve(). 4128 4129 Options Database Key: 4130 . -ts_max_snes_failures - Maximum number of nonlinear solve failures 4131 4132 Level: intermediate 4133 4134 .keywords: TS, set, maximum, number 4135 4136 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason() 4137 @*/ 4138 PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails) 4139 { 4140 PetscFunctionBegin; 4141 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4142 ts->max_snes_failures = fails; 4143 PetscFunctionReturn(0); 4144 } 4145 4146 #undef __FUNCT__ 4147 #define __FUNCT__ "TSSetErrorIfStepFails" 4148 /*@ 4149 TSSetErrorIfStepFails - Error if no step succeeds 4150 4151 Not Collective 4152 4153 Input Parameter: 4154 + ts - TS context 4155 - err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure 4156 4157 Options Database Key: 4158 . -ts_error_if_step_fails - Error if no step succeeds 4159 4160 Level: intermediate 4161 4162 .keywords: TS, set, error 4163 4164 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 4165 @*/ 4166 PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err) 4167 { 4168 PetscFunctionBegin; 4169 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4170 ts->errorifstepfailed = err; 4171 PetscFunctionReturn(0); 4172 } 4173 4174 #undef __FUNCT__ 4175 #define __FUNCT__ "TSMonitorSolutionBinary" 4176 /*@C 4177 TSMonitorSolutionBinary - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file 4178 4179 Collective on TS 4180 4181 Input Parameters: 4182 + ts - the TS context 4183 . step - current time-step 4184 . ptime - current time 4185 . u - current state 4186 - viewer - binary viewer 4187 4188 Level: intermediate 4189 4190 .keywords: TS, vector, monitor, view 4191 4192 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4193 @*/ 4194 PetscErrorCode TSMonitorSolutionBinary(TS ts,PetscInt step,PetscReal ptime,Vec u,void *viewer) 4195 { 4196 PetscErrorCode ierr; 4197 PetscViewer v = (PetscViewer)viewer; 4198 4199 PetscFunctionBegin; 4200 ierr = VecView(u,v);CHKERRQ(ierr); 4201 PetscFunctionReturn(0); 4202 } 4203 4204 #undef __FUNCT__ 4205 #define __FUNCT__ "TSMonitorSolutionVTK" 4206 /*@C 4207 TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep. 4208 4209 Collective on TS 4210 4211 Input Parameters: 4212 + ts - the TS context 4213 . step - current time-step 4214 . ptime - current time 4215 . u - current state 4216 - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 4217 4218 Level: intermediate 4219 4220 Notes: 4221 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. 4222 These are named according to the file name template. 4223 4224 This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy(). 4225 4226 .keywords: TS, vector, monitor, view 4227 4228 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4229 @*/ 4230 PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec u,void *filenametemplate) 4231 { 4232 PetscErrorCode ierr; 4233 char filename[PETSC_MAX_PATH_LEN]; 4234 PetscViewer viewer; 4235 4236 PetscFunctionBegin; 4237 ierr = PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);CHKERRQ(ierr); 4238 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)ts),filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 4239 ierr = VecView(u,viewer);CHKERRQ(ierr); 4240 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 4241 PetscFunctionReturn(0); 4242 } 4243 4244 #undef __FUNCT__ 4245 #define __FUNCT__ "TSMonitorSolutionVTKDestroy" 4246 /*@C 4247 TSMonitorSolutionVTKDestroy - Destroy context for monitoring 4248 4249 Collective on TS 4250 4251 Input Parameters: 4252 . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 4253 4254 Level: intermediate 4255 4256 Note: 4257 This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK(). 4258 4259 .keywords: TS, vector, monitor, view 4260 4261 .seealso: TSMonitorSet(), TSMonitorSolutionVTK() 4262 @*/ 4263 PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate) 4264 { 4265 PetscErrorCode ierr; 4266 4267 PetscFunctionBegin; 4268 ierr = PetscFree(*(char**)filenametemplate);CHKERRQ(ierr); 4269 PetscFunctionReturn(0); 4270 } 4271 4272 #undef __FUNCT__ 4273 #define __FUNCT__ "TSGetAdapt" 4274 /*@ 4275 TSGetAdapt - Get the adaptive controller context for the current method 4276 4277 Collective on TS if controller has not been created yet 4278 4279 Input Arguments: 4280 . ts - time stepping context 4281 4282 Output Arguments: 4283 . adapt - adaptive controller 4284 4285 Level: intermediate 4286 4287 .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose() 4288 @*/ 4289 PetscErrorCode TSGetAdapt(TS ts,TSAdapt *adapt) 4290 { 4291 PetscErrorCode ierr; 4292 4293 PetscFunctionBegin; 4294 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4295 PetscValidPointer(adapt,2); 4296 if (!ts->adapt) { 4297 ierr = TSAdaptCreate(PetscObjectComm((PetscObject)ts),&ts->adapt);CHKERRQ(ierr); 4298 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->adapt);CHKERRQ(ierr); 4299 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);CHKERRQ(ierr); 4300 } 4301 *adapt = ts->adapt; 4302 PetscFunctionReturn(0); 4303 } 4304 4305 #undef __FUNCT__ 4306 #define __FUNCT__ "TSSetTolerances" 4307 /*@ 4308 TSSetTolerances - Set tolerances for local truncation error when using adaptive controller 4309 4310 Logically Collective 4311 4312 Input Arguments: 4313 + ts - time integration context 4314 . atol - scalar absolute tolerances, PETSC_DECIDE to leave current value 4315 . vatol - vector of absolute tolerances or NULL, used in preference to atol if present 4316 . rtol - scalar relative tolerances, PETSC_DECIDE to leave current value 4317 - vrtol - vector of relative tolerances or NULL, used in preference to atol if present 4318 4319 Options Database keys: 4320 + -ts_rtol <rtol> - relative tolerance for local truncation error 4321 - -ts_atol <atol> Absolute tolerance for local truncation error 4322 4323 Level: beginner 4324 4325 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSGetTolerances() 4326 @*/ 4327 PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol) 4328 { 4329 PetscErrorCode ierr; 4330 4331 PetscFunctionBegin; 4332 if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol; 4333 if (vatol) { 4334 ierr = PetscObjectReference((PetscObject)vatol);CHKERRQ(ierr); 4335 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 4336 4337 ts->vatol = vatol; 4338 } 4339 if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol; 4340 if (vrtol) { 4341 ierr = PetscObjectReference((PetscObject)vrtol);CHKERRQ(ierr); 4342 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 4343 4344 ts->vrtol = vrtol; 4345 } 4346 PetscFunctionReturn(0); 4347 } 4348 4349 #undef __FUNCT__ 4350 #define __FUNCT__ "TSGetTolerances" 4351 /*@ 4352 TSGetTolerances - Get tolerances for local truncation error when using adaptive controller 4353 4354 Logically Collective 4355 4356 Input Arguments: 4357 . ts - time integration context 4358 4359 Output Arguments: 4360 + atol - scalar absolute tolerances, NULL to ignore 4361 . vatol - vector of absolute tolerances, NULL to ignore 4362 . rtol - scalar relative tolerances, NULL to ignore 4363 - vrtol - vector of relative tolerances, NULL to ignore 4364 4365 Level: beginner 4366 4367 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSSetTolerances() 4368 @*/ 4369 PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol) 4370 { 4371 PetscFunctionBegin; 4372 if (atol) *atol = ts->atol; 4373 if (vatol) *vatol = ts->vatol; 4374 if (rtol) *rtol = ts->rtol; 4375 if (vrtol) *vrtol = ts->vrtol; 4376 PetscFunctionReturn(0); 4377 } 4378 4379 #undef __FUNCT__ 4380 #define __FUNCT__ "TSErrorNormWRMS" 4381 /*@ 4382 TSErrorNormWRMS - compute a weighted norm of the difference between a vector and the current state 4383 4384 Collective on TS 4385 4386 Input Arguments: 4387 + ts - time stepping context 4388 - Y - state vector to be compared to ts->vec_sol 4389 4390 Output Arguments: 4391 . norm - weighted norm, a value of 1.0 is considered small 4392 4393 Level: developer 4394 4395 .seealso: TSSetTolerances() 4396 @*/ 4397 PetscErrorCode TSErrorNormWRMS(TS ts,Vec Y,PetscReal *norm) 4398 { 4399 PetscErrorCode ierr; 4400 PetscInt i,n,N; 4401 const PetscScalar *u,*y; 4402 Vec U; 4403 PetscReal sum,gsum; 4404 4405 PetscFunctionBegin; 4406 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4407 PetscValidHeaderSpecific(Y,VEC_CLASSID,2); 4408 PetscValidPointer(norm,3); 4409 U = ts->vec_sol; 4410 PetscCheckSameTypeAndComm(U,1,Y,2); 4411 if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"Y cannot be the TS solution vector"); 4412 4413 ierr = VecGetSize(U,&N);CHKERRQ(ierr); 4414 ierr = VecGetLocalSize(U,&n);CHKERRQ(ierr); 4415 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 4416 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 4417 sum = 0.; 4418 if (ts->vatol && ts->vrtol) { 4419 const PetscScalar *atol,*rtol; 4420 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4421 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4422 for (i=0; i<n; i++) { 4423 PetscReal tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4424 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4425 } 4426 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4427 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4428 } else if (ts->vatol) { /* vector atol, scalar rtol */ 4429 const PetscScalar *atol; 4430 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4431 for (i=0; i<n; i++) { 4432 PetscReal tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4433 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4434 } 4435 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4436 } else if (ts->vrtol) { /* scalar atol, vector rtol */ 4437 const PetscScalar *rtol; 4438 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4439 for (i=0; i<n; i++) { 4440 PetscReal tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4441 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4442 } 4443 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4444 } else { /* scalar atol, scalar rtol */ 4445 for (i=0; i<n; i++) { 4446 PetscReal tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4447 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4448 } 4449 } 4450 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 4451 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 4452 4453 ierr = MPI_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 4454 *norm = PetscSqrtReal(gsum / N); 4455 if (PetscIsInfOrNanReal(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 4456 PetscFunctionReturn(0); 4457 } 4458 4459 #undef __FUNCT__ 4460 #define __FUNCT__ "TSSetCFLTimeLocal" 4461 /*@ 4462 TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler 4463 4464 Logically Collective on TS 4465 4466 Input Arguments: 4467 + ts - time stepping context 4468 - cfltime - maximum stable time step if using forward Euler (value can be different on each process) 4469 4470 Note: 4471 After calling this function, the global CFL time can be obtained by calling TSGetCFLTime() 4472 4473 Level: intermediate 4474 4475 .seealso: TSGetCFLTime(), TSADAPTCFL 4476 @*/ 4477 PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime) 4478 { 4479 PetscFunctionBegin; 4480 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4481 ts->cfltime_local = cfltime; 4482 ts->cfltime = -1.; 4483 PetscFunctionReturn(0); 4484 } 4485 4486 #undef __FUNCT__ 4487 #define __FUNCT__ "TSGetCFLTime" 4488 /*@ 4489 TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler 4490 4491 Collective on TS 4492 4493 Input Arguments: 4494 . ts - time stepping context 4495 4496 Output Arguments: 4497 . cfltime - maximum stable time step for forward Euler 4498 4499 Level: advanced 4500 4501 .seealso: TSSetCFLTimeLocal() 4502 @*/ 4503 PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime) 4504 { 4505 PetscErrorCode ierr; 4506 4507 PetscFunctionBegin; 4508 if (ts->cfltime < 0) { 4509 ierr = MPI_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 4510 } 4511 *cfltime = ts->cfltime; 4512 PetscFunctionReturn(0); 4513 } 4514 4515 #undef __FUNCT__ 4516 #define __FUNCT__ "TSVISetVariableBounds" 4517 /*@ 4518 TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu 4519 4520 Input Parameters: 4521 . ts - the TS context. 4522 . xl - lower bound. 4523 . xu - upper bound. 4524 4525 Notes: 4526 If this routine is not called then the lower and upper bounds are set to 4527 PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp(). 4528 4529 Level: advanced 4530 4531 @*/ 4532 PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu) 4533 { 4534 PetscErrorCode ierr; 4535 SNES snes; 4536 4537 PetscFunctionBegin; 4538 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 4539 ierr = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(ierr); 4540 PetscFunctionReturn(0); 4541 } 4542 4543 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4544 #include <mex.h> 4545 4546 typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext; 4547 4548 #undef __FUNCT__ 4549 #define __FUNCT__ "TSComputeFunction_Matlab" 4550 /* 4551 TSComputeFunction_Matlab - Calls the function that has been set with 4552 TSSetFunctionMatlab(). 4553 4554 Collective on TS 4555 4556 Input Parameters: 4557 + snes - the TS context 4558 - u - input vector 4559 4560 Output Parameter: 4561 . y - function vector, as set by TSSetFunction() 4562 4563 Notes: 4564 TSComputeFunction() is typically used within nonlinear solvers 4565 implementations, so most users would not generally call this routine 4566 themselves. 4567 4568 Level: developer 4569 4570 .keywords: TS, nonlinear, compute, function 4571 4572 .seealso: TSSetFunction(), TSGetFunction() 4573 */ 4574 PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec u,Vec udot,Vec y, void *ctx) 4575 { 4576 PetscErrorCode ierr; 4577 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 4578 int nlhs = 1,nrhs = 7; 4579 mxArray *plhs[1],*prhs[7]; 4580 long long int lx = 0,lxdot = 0,ly = 0,ls = 0; 4581 4582 PetscFunctionBegin; 4583 PetscValidHeaderSpecific(snes,TS_CLASSID,1); 4584 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 4585 PetscValidHeaderSpecific(udot,VEC_CLASSID,4); 4586 PetscValidHeaderSpecific(y,VEC_CLASSID,5); 4587 PetscCheckSameComm(snes,1,u,3); 4588 PetscCheckSameComm(snes,1,y,5); 4589 4590 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4591 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4592 ierr = PetscMemcpy(&lxdot,&udot,sizeof(udot));CHKERRQ(ierr); 4593 ierr = PetscMemcpy(&ly,&y,sizeof(u));CHKERRQ(ierr); 4594 4595 prhs[0] = mxCreateDoubleScalar((double)ls); 4596 prhs[1] = mxCreateDoubleScalar(time); 4597 prhs[2] = mxCreateDoubleScalar((double)lx); 4598 prhs[3] = mxCreateDoubleScalar((double)lxdot); 4599 prhs[4] = mxCreateDoubleScalar((double)ly); 4600 prhs[5] = mxCreateString(sctx->funcname); 4601 prhs[6] = sctx->ctx; 4602 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");CHKERRQ(ierr); 4603 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4604 mxDestroyArray(prhs[0]); 4605 mxDestroyArray(prhs[1]); 4606 mxDestroyArray(prhs[2]); 4607 mxDestroyArray(prhs[3]); 4608 mxDestroyArray(prhs[4]); 4609 mxDestroyArray(prhs[5]); 4610 mxDestroyArray(plhs[0]); 4611 PetscFunctionReturn(0); 4612 } 4613 4614 4615 #undef __FUNCT__ 4616 #define __FUNCT__ "TSSetFunctionMatlab" 4617 /* 4618 TSSetFunctionMatlab - Sets the function evaluation routine and function 4619 vector for use by the TS routines in solving ODEs 4620 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4621 4622 Logically Collective on TS 4623 4624 Input Parameters: 4625 + ts - the TS context 4626 - func - function evaluation routine 4627 4628 Calling sequence of func: 4629 $ func (TS ts,PetscReal time,Vec u,Vec udot,Vec f,void *ctx); 4630 4631 Level: beginner 4632 4633 .keywords: TS, nonlinear, set, function 4634 4635 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4636 */ 4637 PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx) 4638 { 4639 PetscErrorCode ierr; 4640 TSMatlabContext *sctx; 4641 4642 PetscFunctionBegin; 4643 /* currently sctx is memory bleed */ 4644 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4645 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4646 /* 4647 This should work, but it doesn't 4648 sctx->ctx = ctx; 4649 mexMakeArrayPersistent(sctx->ctx); 4650 */ 4651 sctx->ctx = mxDuplicateArray(ctx); 4652 4653 ierr = TSSetIFunction(ts,NULL,TSComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4654 PetscFunctionReturn(0); 4655 } 4656 4657 #undef __FUNCT__ 4658 #define __FUNCT__ "TSComputeJacobian_Matlab" 4659 /* 4660 TSComputeJacobian_Matlab - Calls the function that has been set with 4661 TSSetJacobianMatlab(). 4662 4663 Collective on TS 4664 4665 Input Parameters: 4666 + ts - the TS context 4667 . u - input vector 4668 . A, B - the matrices 4669 - ctx - user context 4670 4671 Level: developer 4672 4673 .keywords: TS, nonlinear, compute, function 4674 4675 .seealso: TSSetFunction(), TSGetFunction() 4676 @*/ 4677 PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec u,Vec udot,PetscReal shift,Mat A,Mat B,void *ctx) 4678 { 4679 PetscErrorCode ierr; 4680 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 4681 int nlhs = 2,nrhs = 9; 4682 mxArray *plhs[2],*prhs[9]; 4683 long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0; 4684 4685 PetscFunctionBegin; 4686 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4687 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 4688 4689 /* call Matlab function in ctx with arguments u and y */ 4690 4691 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 4692 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4693 ierr = PetscMemcpy(&lxdot,&udot,sizeof(u));CHKERRQ(ierr); 4694 ierr = PetscMemcpy(&lA,A,sizeof(u));CHKERRQ(ierr); 4695 ierr = PetscMemcpy(&lB,B,sizeof(u));CHKERRQ(ierr); 4696 4697 prhs[0] = mxCreateDoubleScalar((double)ls); 4698 prhs[1] = mxCreateDoubleScalar((double)time); 4699 prhs[2] = mxCreateDoubleScalar((double)lx); 4700 prhs[3] = mxCreateDoubleScalar((double)lxdot); 4701 prhs[4] = mxCreateDoubleScalar((double)shift); 4702 prhs[5] = mxCreateDoubleScalar((double)lA); 4703 prhs[6] = mxCreateDoubleScalar((double)lB); 4704 prhs[7] = mxCreateString(sctx->funcname); 4705 prhs[8] = sctx->ctx; 4706 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");CHKERRQ(ierr); 4707 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4708 mxDestroyArray(prhs[0]); 4709 mxDestroyArray(prhs[1]); 4710 mxDestroyArray(prhs[2]); 4711 mxDestroyArray(prhs[3]); 4712 mxDestroyArray(prhs[4]); 4713 mxDestroyArray(prhs[5]); 4714 mxDestroyArray(prhs[6]); 4715 mxDestroyArray(prhs[7]); 4716 mxDestroyArray(plhs[0]); 4717 mxDestroyArray(plhs[1]); 4718 PetscFunctionReturn(0); 4719 } 4720 4721 4722 #undef __FUNCT__ 4723 #define __FUNCT__ "TSSetJacobianMatlab" 4724 /* 4725 TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4726 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 4727 4728 Logically Collective on TS 4729 4730 Input Parameters: 4731 + ts - the TS context 4732 . A,B - Jacobian matrices 4733 . func - function evaluation routine 4734 - ctx - user context 4735 4736 Calling sequence of func: 4737 $ flag = func (TS ts,PetscReal time,Vec u,Vec udot,Mat A,Mat B,void *ctx); 4738 4739 4740 Level: developer 4741 4742 .keywords: TS, nonlinear, set, function 4743 4744 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4745 */ 4746 PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx) 4747 { 4748 PetscErrorCode ierr; 4749 TSMatlabContext *sctx; 4750 4751 PetscFunctionBegin; 4752 /* currently sctx is memory bleed */ 4753 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4754 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4755 /* 4756 This should work, but it doesn't 4757 sctx->ctx = ctx; 4758 mexMakeArrayPersistent(sctx->ctx); 4759 */ 4760 sctx->ctx = mxDuplicateArray(ctx); 4761 4762 ierr = TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4763 PetscFunctionReturn(0); 4764 } 4765 4766 #undef __FUNCT__ 4767 #define __FUNCT__ "TSMonitor_Matlab" 4768 /* 4769 TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab(). 4770 4771 Collective on TS 4772 4773 .seealso: TSSetFunction(), TSGetFunction() 4774 @*/ 4775 PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec u, void *ctx) 4776 { 4777 PetscErrorCode ierr; 4778 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 4779 int nlhs = 1,nrhs = 6; 4780 mxArray *plhs[1],*prhs[6]; 4781 long long int lx = 0,ls = 0; 4782 4783 PetscFunctionBegin; 4784 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4785 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 4786 4787 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 4788 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4789 4790 prhs[0] = mxCreateDoubleScalar((double)ls); 4791 prhs[1] = mxCreateDoubleScalar((double)it); 4792 prhs[2] = mxCreateDoubleScalar((double)time); 4793 prhs[3] = mxCreateDoubleScalar((double)lx); 4794 prhs[4] = mxCreateString(sctx->funcname); 4795 prhs[5] = sctx->ctx; 4796 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");CHKERRQ(ierr); 4797 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4798 mxDestroyArray(prhs[0]); 4799 mxDestroyArray(prhs[1]); 4800 mxDestroyArray(prhs[2]); 4801 mxDestroyArray(prhs[3]); 4802 mxDestroyArray(prhs[4]); 4803 mxDestroyArray(plhs[0]); 4804 PetscFunctionReturn(0); 4805 } 4806 4807 4808 #undef __FUNCT__ 4809 #define __FUNCT__ "TSMonitorSetMatlab" 4810 /* 4811 TSMonitorSetMatlab - Sets the monitor function from Matlab 4812 4813 Level: developer 4814 4815 .keywords: TS, nonlinear, set, function 4816 4817 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4818 */ 4819 PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx) 4820 { 4821 PetscErrorCode ierr; 4822 TSMatlabContext *sctx; 4823 4824 PetscFunctionBegin; 4825 /* currently sctx is memory bleed */ 4826 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4827 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4828 /* 4829 This should work, but it doesn't 4830 sctx->ctx = ctx; 4831 mexMakeArrayPersistent(sctx->ctx); 4832 */ 4833 sctx->ctx = mxDuplicateArray(ctx); 4834 4835 ierr = TSMonitorSet(ts,TSMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 4836 PetscFunctionReturn(0); 4837 } 4838 #endif 4839 4840 4841 4842 #undef __FUNCT__ 4843 #define __FUNCT__ "TSMonitorLGSolution" 4844 /*@C 4845 TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector 4846 in a time based line graph 4847 4848 Collective on TS 4849 4850 Input Parameters: 4851 + ts - the TS context 4852 . step - current time-step 4853 . ptime - current time 4854 - lg - a line graph object 4855 4856 Level: intermediate 4857 4858 Notes: each process in a parallel run displays its component solutions in a separate window 4859 4860 .keywords: TS, vector, monitor, view 4861 4862 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4863 @*/ 4864 PetscErrorCode TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4865 { 4866 PetscErrorCode ierr; 4867 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 4868 const PetscScalar *yy; 4869 PetscInt dim; 4870 4871 PetscFunctionBegin; 4872 if (!step) { 4873 PetscDrawAxis axis; 4874 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4875 ierr = PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");CHKERRQ(ierr); 4876 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 4877 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 4878 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4879 } 4880 ierr = VecGetArrayRead(u,&yy);CHKERRQ(ierr); 4881 #if defined(PETSC_USE_COMPLEX) 4882 { 4883 PetscReal *yreal; 4884 PetscInt i,n; 4885 ierr = VecGetLocalSize(u,&n);CHKERRQ(ierr); 4886 ierr = PetscMalloc1(n,&yreal);CHKERRQ(ierr); 4887 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 4888 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 4889 ierr = PetscFree(yreal);CHKERRQ(ierr); 4890 } 4891 #else 4892 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 4893 #endif 4894 ierr = VecRestoreArrayRead(u,&yy);CHKERRQ(ierr); 4895 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 4896 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4897 } 4898 PetscFunctionReturn(0); 4899 } 4900 4901 #undef __FUNCT__ 4902 #define __FUNCT__ "TSMonitorLGError" 4903 /*@C 4904 TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the solution vector 4905 in a time based line graph 4906 4907 Collective on TS 4908 4909 Input Parameters: 4910 + ts - the TS context 4911 . step - current time-step 4912 . ptime - current time 4913 - lg - a line graph object 4914 4915 Level: intermediate 4916 4917 Notes: 4918 Only for sequential solves. 4919 4920 The user must provide the solution using TSSetSolutionFunction() to use this monitor. 4921 4922 Options Database Keys: 4923 . -ts_monitor_lg_error - create a graphical monitor of error history 4924 4925 .keywords: TS, vector, monitor, view 4926 4927 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction() 4928 @*/ 4929 PetscErrorCode TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4930 { 4931 PetscErrorCode ierr; 4932 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 4933 const PetscScalar *yy; 4934 Vec y; 4935 PetscInt dim; 4936 4937 PetscFunctionBegin; 4938 if (!step) { 4939 PetscDrawAxis axis; 4940 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4941 ierr = PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Solution");CHKERRQ(ierr); 4942 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 4943 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 4944 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4945 } 4946 ierr = VecDuplicate(u,&y);CHKERRQ(ierr); 4947 ierr = TSComputeSolutionFunction(ts,ptime,y);CHKERRQ(ierr); 4948 ierr = VecAXPY(y,-1.0,u);CHKERRQ(ierr); 4949 ierr = VecGetArrayRead(y,&yy);CHKERRQ(ierr); 4950 #if defined(PETSC_USE_COMPLEX) 4951 { 4952 PetscReal *yreal; 4953 PetscInt i,n; 4954 ierr = VecGetLocalSize(y,&n);CHKERRQ(ierr); 4955 ierr = PetscMalloc1(n,&yreal);CHKERRQ(ierr); 4956 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 4957 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 4958 ierr = PetscFree(yreal);CHKERRQ(ierr); 4959 } 4960 #else 4961 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 4962 #endif 4963 ierr = VecRestoreArrayRead(y,&yy);CHKERRQ(ierr); 4964 ierr = VecDestroy(&y);CHKERRQ(ierr); 4965 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 4966 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4967 } 4968 PetscFunctionReturn(0); 4969 } 4970 4971 #undef __FUNCT__ 4972 #define __FUNCT__ "TSMonitorLGSNESIterations" 4973 PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 4974 { 4975 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 4976 PetscReal x = ptime,y; 4977 PetscErrorCode ierr; 4978 PetscInt its; 4979 4980 PetscFunctionBegin; 4981 if (!n) { 4982 PetscDrawAxis axis; 4983 4984 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4985 ierr = PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");CHKERRQ(ierr); 4986 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4987 4988 ctx->snes_its = 0; 4989 } 4990 ierr = TSGetSNESIterations(ts,&its);CHKERRQ(ierr); 4991 y = its - ctx->snes_its; 4992 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 4993 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 4994 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4995 } 4996 ctx->snes_its = its; 4997 PetscFunctionReturn(0); 4998 } 4999 5000 #undef __FUNCT__ 5001 #define __FUNCT__ "TSMonitorLGKSPIterations" 5002 PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 5003 { 5004 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 5005 PetscReal x = ptime,y; 5006 PetscErrorCode ierr; 5007 PetscInt its; 5008 5009 PetscFunctionBegin; 5010 if (!n) { 5011 PetscDrawAxis axis; 5012 5013 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 5014 ierr = PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");CHKERRQ(ierr); 5015 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 5016 5017 ctx->ksp_its = 0; 5018 } 5019 ierr = TSGetKSPIterations(ts,&its);CHKERRQ(ierr); 5020 y = its - ctx->ksp_its; 5021 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 5022 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 5023 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 5024 } 5025 ctx->ksp_its = its; 5026 PetscFunctionReturn(0); 5027 } 5028 5029 #undef __FUNCT__ 5030 #define __FUNCT__ "TSComputeLinearStability" 5031 /*@ 5032 TSComputeLinearStability - computes the linear stability function at a point 5033 5034 Collective on TS and Vec 5035 5036 Input Parameters: 5037 + ts - the TS context 5038 - xr,xi - real and imaginary part of input arguments 5039 5040 Output Parameters: 5041 . yr,yi - real and imaginary part of function value 5042 5043 Level: developer 5044 5045 .keywords: TS, compute 5046 5047 .seealso: TSSetRHSFunction(), TSComputeIFunction() 5048 @*/ 5049 PetscErrorCode TSComputeLinearStability(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi) 5050 { 5051 PetscErrorCode ierr; 5052 5053 PetscFunctionBegin; 5054 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5055 if (!ts->ops->linearstability) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Linearized stability function not provided for this method"); 5056 ierr = (*ts->ops->linearstability)(ts,xr,xi,yr,yi);CHKERRQ(ierr); 5057 PetscFunctionReturn(0); 5058 } 5059 5060 #undef __FUNCT__ 5061 #define __FUNCT__ "TSRollBack" 5062 /*@ 5063 TSRollBack - Rolls back one time step 5064 5065 Collective on TS 5066 5067 Input Parameter: 5068 . ts - the TS context obtained from TSCreate() 5069 5070 Level: advanced 5071 5072 .keywords: TS, timestep, rollback 5073 5074 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSInterpolate() 5075 @*/ 5076 PetscErrorCode TSRollBack(TS ts) 5077 { 5078 PetscErrorCode ierr; 5079 5080 PetscFunctionBegin; 5081 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 5082 5083 if (!ts->ops->rollback) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRollBack not implemented for type '%s'",((PetscObject)ts)->type_name); 5084 ierr = (*ts->ops->rollback)(ts);CHKERRQ(ierr); 5085 ts->time_step = ts->ptime - ts->ptime_prev; 5086 ts->ptime = ts->ptime_prev; 5087 PetscFunctionReturn(0); 5088 } 5089 5090 #undef __FUNCT__ 5091 #define __FUNCT__ "TSGetStages" 5092 /*@ 5093 TSGetStages - Get the number of stages and stage values 5094 5095 Input Parameter: 5096 . ts - the TS context obtained from TSCreate() 5097 5098 Level: advanced 5099 5100 .keywords: TS, getstages 5101 5102 .seealso: TSCreate() 5103 @*/ 5104 PetscErrorCode TSGetStages(TS ts,PetscInt *ns, Vec **Y) 5105 { 5106 PetscErrorCode ierr; 5107 5108 PetscFunctionBegin; 5109 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 5110 PetscValidPointer(ns,2); 5111 5112 if (!ts->ops->getstages) *ns=0; 5113 else { 5114 ierr = (*ts->ops->getstages)(ts,ns,Y);CHKERRQ(ierr); 5115 } 5116 PetscFunctionReturn(0); 5117 } 5118 5119 5120 5121