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