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