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