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