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