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