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