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