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_use_markers <true,false> - mark 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 = PetscDrawLGSetUseMarkers((*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 = PetscDrawLGSetUseMarkers(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__ "TSErrorWeightedNorm2" 4967 /*@ 4968 TSErrorWeightedNorm2 - compute a weighted 2-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 Output Arguments: 4977 . norm - weighted norm, a value of 1.0 is considered small 4978 4979 Level: developer 4980 4981 .seealso: TSErrorWeightedNorm(), TSErrorWeightedNormInfinity() 4982 @*/ 4983 PetscErrorCode TSErrorWeightedNorm2(TS ts,Vec Y,PetscReal *norm) 4984 { 4985 PetscErrorCode ierr; 4986 PetscInt i,n,N,rstart; 4987 const PetscScalar *u,*y; 4988 Vec U; 4989 PetscReal sum,gsum; 4990 PetscReal tol; 4991 4992 PetscFunctionBegin; 4993 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4994 PetscValidHeaderSpecific(Y,VEC_CLASSID,2); 4995 PetscValidPointer(norm,3); 4996 U = ts->vec_sol; 4997 PetscCheckSameTypeAndComm(U,1,Y,2); 4998 if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"Y cannot be the TS solution vector"); 4999 5000 ierr = VecGetSize(U,&N);CHKERRQ(ierr); 5001 ierr = VecGetLocalSize(U,&n);CHKERRQ(ierr); 5002 ierr = VecGetOwnershipRange(U,&rstart,NULL);CHKERRQ(ierr); 5003 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 5004 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 5005 sum = 0.; 5006 if (ts->vatol && ts->vrtol) { 5007 const PetscScalar *atol,*rtol; 5008 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 5009 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 5010 if(ts->is_diff) { 5011 const PetscInt *idx; 5012 PetscInt k; 5013 ierr = ISGetIndices(ts->is_diff,&idx);CHKERRQ(ierr); 5014 ierr = ISGetLocalSize(ts->is_diff,&n);CHKERRQ(ierr); 5015 for(i=0; i < n; i++) { 5016 k = idx[i] - rstart; 5017 tol = PetscRealPart(atol[k]) + PetscRealPart(rtol[k])*PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5018 sum += PetscSqr(PetscAbsScalar(y[k] - u[k]) / tol); 5019 } 5020 ierr = ISRestoreIndices(ts->is_diff,&idx);CHKERRQ(ierr); 5021 } else { 5022 for (i=0; i<n; i++) { 5023 tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 5024 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 5025 } 5026 } 5027 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 5028 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 5029 } else if (ts->vatol) { /* vector atol, scalar rtol */ 5030 const PetscScalar *atol; 5031 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 5032 if(ts->is_diff) { 5033 const PetscInt *idx; 5034 PetscInt k; 5035 ierr = ISGetIndices(ts->is_diff,&idx);CHKERRQ(ierr); 5036 ierr = ISGetLocalSize(ts->is_diff,&n);CHKERRQ(ierr); 5037 for(i=0; i < n; i++) { 5038 k = idx[i] - rstart; 5039 tol = PetscRealPart(atol[k]) + ts->rtol * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5040 sum += PetscSqr(PetscAbsScalar(y[k] - u[k]) / tol); 5041 } 5042 ierr = ISRestoreIndices(ts->is_diff,&idx);CHKERRQ(ierr); 5043 } else { 5044 for (i=0; i<n; i++) { 5045 tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 5046 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 5047 } 5048 } 5049 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 5050 } else if (ts->vrtol) { /* scalar atol, vector rtol */ 5051 const PetscScalar *rtol; 5052 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 5053 if(ts->is_diff) { 5054 const PetscInt *idx; 5055 PetscInt k; 5056 ierr = ISGetIndices(ts->is_diff,&idx);CHKERRQ(ierr); 5057 ierr = ISGetLocalSize(ts->is_diff,&n);CHKERRQ(ierr); 5058 for(i=0; i < n; i++) { 5059 k = idx[i] - rstart; 5060 tol = ts->atol + PetscRealPart(rtol[k]) * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5061 sum += PetscSqr(PetscAbsScalar(y[k] - u[k]) / tol); 5062 } 5063 ierr = ISRestoreIndices(ts->is_diff,&idx);CHKERRQ(ierr); 5064 } else { 5065 for (i=0; i<n; i++) { 5066 tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 5067 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 5068 } 5069 } 5070 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 5071 } else { /* scalar atol, scalar rtol */ 5072 if (ts->is_diff) { 5073 const PetscInt *idx; 5074 PetscInt k; 5075 ierr = ISGetIndices(ts->is_diff,&idx);CHKERRQ(ierr); 5076 ierr = ISGetLocalSize(ts->is_diff,&n);CHKERRQ(ierr); 5077 for (i=0; i<n; i++) { 5078 k = idx[i] - rstart; 5079 tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5080 sum += PetscSqr(PetscAbsScalar(y[k] - u[k]) / tol); 5081 } 5082 } else { 5083 for (i=0; i<n; i++) { 5084 tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 5085 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 5086 } 5087 } 5088 } 5089 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 5090 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 5091 5092 ierr = MPI_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 5093 *norm = PetscSqrtReal(gsum / N); 5094 5095 if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 5096 PetscFunctionReturn(0); 5097 } 5098 5099 #undef __FUNCT__ 5100 #define __FUNCT__ "TSErrorWeightedNormInfinity" 5101 /*@ 5102 TSErrorWeightedNormInfinity - compute a weighted infinity-norm of the difference between a vector and the current state 5103 5104 Collective on TS 5105 5106 Input Arguments: 5107 + ts - time stepping context 5108 - Y - state vector to be compared to ts->vec_sol 5109 5110 Output Arguments: 5111 . norm - weighted norm, a value of 1.0 is considered small 5112 5113 Level: developer 5114 5115 .seealso: TSErrorWeightedNorm(), TSErrorWeightedNorm2() 5116 @*/ 5117 PetscErrorCode TSErrorWeightedNormInfinity(TS ts,Vec Y,PetscReal *norm) 5118 { 5119 PetscErrorCode ierr; 5120 PetscInt i,n,N,rstart,k; 5121 const PetscScalar *u,*y; 5122 Vec U; 5123 PetscReal max,gmax; 5124 PetscReal tol; 5125 5126 PetscFunctionBegin; 5127 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5128 PetscValidHeaderSpecific(Y,VEC_CLASSID,2); 5129 PetscValidPointer(norm,3); 5130 U = ts->vec_sol; 5131 PetscCheckSameTypeAndComm(U,1,Y,2); 5132 if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"Y cannot be the TS solution vector"); 5133 5134 ierr = VecGetSize(U,&N);CHKERRQ(ierr); 5135 ierr = VecGetLocalSize(U,&n);CHKERRQ(ierr); 5136 ierr = VecGetOwnershipRange(U,&rstart,NULL);CHKERRQ(ierr); 5137 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 5138 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 5139 if (ts->vatol && ts->vrtol) { 5140 const PetscScalar *atol,*rtol; 5141 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 5142 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 5143 if(ts->is_diff) { 5144 const PetscInt *idx; 5145 ierr = ISGetIndices(ts->is_diff,&idx);CHKERRQ(ierr); 5146 ierr = ISGetLocalSize(ts->is_diff,&n);CHKERRQ(ierr); 5147 5148 k = idx[0]; 5149 tol = PetscRealPart(atol[k]) + PetscRealPart(rtol[k]) * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5150 max = PetscAbsScalar(y[k] - u[k]) / tol; 5151 for(i=1; i < n; i++) { 5152 k = idx[i] - rstart; 5153 tol = PetscRealPart(atol[k]) + PetscRealPart(rtol[k])*PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5154 max = PetscMax(max,PetscAbsScalar(y[k] - u[k]) / tol); 5155 } 5156 ierr = ISRestoreIndices(ts->is_diff,&idx);CHKERRQ(ierr); 5157 } else { 5158 k = 0; 5159 tol = PetscRealPart(atol[k]) + PetscRealPart(rtol[k]) * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5160 max = PetscAbsScalar(y[k] - u[k]) / tol; 5161 for (i=1; i<n; i++) { 5162 tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 5163 max = PetscMax(max,PetscAbsScalar(y[i] - u[i]) / tol); 5164 } 5165 } 5166 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 5167 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 5168 } else if (ts->vatol) { /* vector atol, scalar rtol */ 5169 const PetscScalar *atol; 5170 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 5171 if(ts->is_diff) { 5172 const PetscInt *idx; 5173 ierr = ISGetIndices(ts->is_diff,&idx);CHKERRQ(ierr); 5174 ierr = ISGetLocalSize(ts->is_diff,&n);CHKERRQ(ierr); 5175 5176 k = idx[0]; 5177 tol = PetscRealPart(atol[k]) + ts->rtol * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5178 max = PetscAbsScalar(y[k] - u[k]) / tol; 5179 for(i=1; i < n; i++) { 5180 k = idx[i] - rstart; 5181 tol = PetscRealPart(atol[k]) + ts->rtol * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5182 max = PetscMax(max,PetscAbsScalar(y[k] - u[k]) / tol); 5183 } 5184 ierr = ISRestoreIndices(ts->is_diff,&idx);CHKERRQ(ierr); 5185 } else { 5186 k = 0; 5187 tol = PetscRealPart(atol[k]) + ts->rtol * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5188 max = PetscAbsScalar(y[k] - u[k]) / tol; 5189 for (i=1; i<n; i++) { 5190 tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 5191 max = PetscMax(max,PetscAbsScalar(y[i] - u[i]) / tol); 5192 } 5193 } 5194 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 5195 } else if (ts->vrtol) { /* scalar atol, vector rtol */ 5196 const PetscScalar *rtol; 5197 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 5198 if(ts->is_diff) { 5199 const PetscInt *idx; 5200 ierr = ISGetIndices(ts->is_diff,&idx);CHKERRQ(ierr); 5201 ierr = ISGetLocalSize(ts->is_diff,&n);CHKERRQ(ierr); 5202 5203 k = idx[0]; 5204 tol = ts->atol + PetscRealPart(rtol[k])*PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5205 max = PetscAbsScalar(y[k] - u[k]) / tol; 5206 for(i=1; i < n; i++) { 5207 k = idx[i] - rstart; 5208 tol = ts->atol + PetscRealPart(rtol[k]) * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5209 max = PetscMax(max,PetscAbsScalar(y[k] - u[k]) / tol); 5210 } 5211 ierr = ISRestoreIndices(ts->is_diff,&idx);CHKERRQ(ierr); 5212 } else { 5213 k = 0; 5214 tol = ts->atol + PetscRealPart(rtol[k]) * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5215 max = PetscAbsScalar(y[k] - u[k]) / tol; 5216 for (i=1; i<n; i++) { 5217 tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 5218 max = PetscMax(max,PetscAbsScalar(y[i] - u[i]) / tol); 5219 } 5220 } 5221 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 5222 } else { /* scalar atol, scalar rtol */ 5223 if (ts->is_diff) { 5224 const PetscInt *idx; 5225 ierr = ISGetIndices(ts->is_diff,&idx);CHKERRQ(ierr); 5226 ierr = ISGetLocalSize(ts->is_diff,&n);CHKERRQ(ierr); 5227 5228 k = idx[0]; 5229 tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5230 max = PetscAbsScalar(y[k] - u[k]) / tol; 5231 for (i=1; i<n; i++) { 5232 k = idx[i] - rstart; 5233 tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5234 max = PetscMax(max,PetscAbsScalar(y[k] - u[k]) / tol); 5235 } 5236 ierr = ISRestoreIndices(ts->is_diff,&idx);CHKERRQ(ierr); 5237 } else { 5238 k = 0; 5239 tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5240 max = PetscAbsScalar(y[k] - u[k]) / tol; 5241 for (i=1; i<n; i++) { 5242 tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 5243 max = PetscMax(max,PetscAbsScalar(y[i] - u[i]) / tol); 5244 } 5245 } 5246 } 5247 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 5248 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 5249 5250 ierr = MPI_Allreduce(&max,&gmax,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 5251 *norm = gmax; 5252 5253 if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 5254 PetscFunctionReturn(0); 5255 } 5256 5257 #undef __FUNCT__ 5258 #define __FUNCT__ "TSErrorWeightedNorm" 5259 /*@ 5260 TSErrorWeightedNorm - compute a weighted norm of the difference between a vector and the current state 5261 5262 Collective on TS 5263 5264 Input Arguments: 5265 + ts - time stepping context 5266 - Y - state vector to be compared to ts->vec_sol 5267 5268 Options Database Keys: 5269 . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY 5270 5271 Output Arguments: 5272 . norm - weighted norm, a value of 1.0 is considered small 5273 5274 Level: developer 5275 5276 .seealso: TSErrorWeightedNormInfinity(), TSErrorWeightedNorm2() 5277 @*/ 5278 PetscErrorCode TSErrorWeightedNorm(TS ts,Vec Y,PetscReal *norm) 5279 { 5280 5281 PetscFunctionBegin; 5282 5283 if(ts->adapt->wnormtype == NORM_2) { 5284 PetscErrorCode ierr; 5285 ierr = TSErrorWeightedNorm2(ts,Y,norm); 5286 } else if(ts->adapt->wnormtype == NORM_INFINITY) { 5287 PetscErrorCode ierr; 5288 ierr = TSErrorWeightedNormInfinity(ts,Y,norm); 5289 } 5290 5291 PetscFunctionReturn(0); 5292 } 5293 5294 5295 #undef __FUNCT__ 5296 #define __FUNCT__ "TSSetCFLTimeLocal" 5297 /*@ 5298 TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler 5299 5300 Logically Collective on TS 5301 5302 Input Arguments: 5303 + ts - time stepping context 5304 - cfltime - maximum stable time step if using forward Euler (value can be different on each process) 5305 5306 Note: 5307 After calling this function, the global CFL time can be obtained by calling TSGetCFLTime() 5308 5309 Level: intermediate 5310 5311 .seealso: TSGetCFLTime(), TSADAPTCFL 5312 @*/ 5313 PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime) 5314 { 5315 PetscFunctionBegin; 5316 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5317 ts->cfltime_local = cfltime; 5318 ts->cfltime = -1.; 5319 PetscFunctionReturn(0); 5320 } 5321 5322 #undef __FUNCT__ 5323 #define __FUNCT__ "TSGetCFLTime" 5324 /*@ 5325 TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler 5326 5327 Collective on TS 5328 5329 Input Arguments: 5330 . ts - time stepping context 5331 5332 Output Arguments: 5333 . cfltime - maximum stable time step for forward Euler 5334 5335 Level: advanced 5336 5337 .seealso: TSSetCFLTimeLocal() 5338 @*/ 5339 PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime) 5340 { 5341 PetscErrorCode ierr; 5342 5343 PetscFunctionBegin; 5344 if (ts->cfltime < 0) { 5345 ierr = MPI_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 5346 } 5347 *cfltime = ts->cfltime; 5348 PetscFunctionReturn(0); 5349 } 5350 5351 #undef __FUNCT__ 5352 #define __FUNCT__ "TSVISetVariableBounds" 5353 /*@ 5354 TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu 5355 5356 Input Parameters: 5357 . ts - the TS context. 5358 . xl - lower bound. 5359 . xu - upper bound. 5360 5361 Notes: 5362 If this routine is not called then the lower and upper bounds are set to 5363 PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp(). 5364 5365 Level: advanced 5366 5367 @*/ 5368 PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu) 5369 { 5370 PetscErrorCode ierr; 5371 SNES snes; 5372 5373 PetscFunctionBegin; 5374 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 5375 ierr = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(ierr); 5376 PetscFunctionReturn(0); 5377 } 5378 5379 #if defined(PETSC_HAVE_MATLAB_ENGINE) 5380 #include <mex.h> 5381 5382 typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext; 5383 5384 #undef __FUNCT__ 5385 #define __FUNCT__ "TSComputeFunction_Matlab" 5386 /* 5387 TSComputeFunction_Matlab - Calls the function that has been set with 5388 TSSetFunctionMatlab(). 5389 5390 Collective on TS 5391 5392 Input Parameters: 5393 + snes - the TS context 5394 - u - input vector 5395 5396 Output Parameter: 5397 . y - function vector, as set by TSSetFunction() 5398 5399 Notes: 5400 TSComputeFunction() is typically used within nonlinear solvers 5401 implementations, so most users would not generally call this routine 5402 themselves. 5403 5404 Level: developer 5405 5406 .keywords: TS, nonlinear, compute, function 5407 5408 .seealso: TSSetFunction(), TSGetFunction() 5409 */ 5410 PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec u,Vec udot,Vec y, void *ctx) 5411 { 5412 PetscErrorCode ierr; 5413 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 5414 int nlhs = 1,nrhs = 7; 5415 mxArray *plhs[1],*prhs[7]; 5416 long long int lx = 0,lxdot = 0,ly = 0,ls = 0; 5417 5418 PetscFunctionBegin; 5419 PetscValidHeaderSpecific(snes,TS_CLASSID,1); 5420 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 5421 PetscValidHeaderSpecific(udot,VEC_CLASSID,4); 5422 PetscValidHeaderSpecific(y,VEC_CLASSID,5); 5423 PetscCheckSameComm(snes,1,u,3); 5424 PetscCheckSameComm(snes,1,y,5); 5425 5426 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5427 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 5428 ierr = PetscMemcpy(&lxdot,&udot,sizeof(udot));CHKERRQ(ierr); 5429 ierr = PetscMemcpy(&ly,&y,sizeof(u));CHKERRQ(ierr); 5430 5431 prhs[0] = mxCreateDoubleScalar((double)ls); 5432 prhs[1] = mxCreateDoubleScalar(time); 5433 prhs[2] = mxCreateDoubleScalar((double)lx); 5434 prhs[3] = mxCreateDoubleScalar((double)lxdot); 5435 prhs[4] = mxCreateDoubleScalar((double)ly); 5436 prhs[5] = mxCreateString(sctx->funcname); 5437 prhs[6] = sctx->ctx; 5438 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");CHKERRQ(ierr); 5439 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5440 mxDestroyArray(prhs[0]); 5441 mxDestroyArray(prhs[1]); 5442 mxDestroyArray(prhs[2]); 5443 mxDestroyArray(prhs[3]); 5444 mxDestroyArray(prhs[4]); 5445 mxDestroyArray(prhs[5]); 5446 mxDestroyArray(plhs[0]); 5447 PetscFunctionReturn(0); 5448 } 5449 5450 5451 #undef __FUNCT__ 5452 #define __FUNCT__ "TSSetFunctionMatlab" 5453 /* 5454 TSSetFunctionMatlab - Sets the function evaluation routine and function 5455 vector for use by the TS routines in solving ODEs 5456 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 5457 5458 Logically Collective on TS 5459 5460 Input Parameters: 5461 + ts - the TS context 5462 - func - function evaluation routine 5463 5464 Calling sequence of func: 5465 $ func (TS ts,PetscReal time,Vec u,Vec udot,Vec f,void *ctx); 5466 5467 Level: beginner 5468 5469 .keywords: TS, nonlinear, set, function 5470 5471 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 5472 */ 5473 PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx) 5474 { 5475 PetscErrorCode ierr; 5476 TSMatlabContext *sctx; 5477 5478 PetscFunctionBegin; 5479 /* currently sctx is memory bleed */ 5480 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 5481 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5482 /* 5483 This should work, but it doesn't 5484 sctx->ctx = ctx; 5485 mexMakeArrayPersistent(sctx->ctx); 5486 */ 5487 sctx->ctx = mxDuplicateArray(ctx); 5488 5489 ierr = TSSetIFunction(ts,NULL,TSComputeFunction_Matlab,sctx);CHKERRQ(ierr); 5490 PetscFunctionReturn(0); 5491 } 5492 5493 #undef __FUNCT__ 5494 #define __FUNCT__ "TSComputeJacobian_Matlab" 5495 /* 5496 TSComputeJacobian_Matlab - Calls the function that has been set with 5497 TSSetJacobianMatlab(). 5498 5499 Collective on TS 5500 5501 Input Parameters: 5502 + ts - the TS context 5503 . u - input vector 5504 . A, B - the matrices 5505 - ctx - user context 5506 5507 Level: developer 5508 5509 .keywords: TS, nonlinear, compute, function 5510 5511 .seealso: TSSetFunction(), TSGetFunction() 5512 @*/ 5513 PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec u,Vec udot,PetscReal shift,Mat A,Mat B,void *ctx) 5514 { 5515 PetscErrorCode ierr; 5516 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 5517 int nlhs = 2,nrhs = 9; 5518 mxArray *plhs[2],*prhs[9]; 5519 long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0; 5520 5521 PetscFunctionBegin; 5522 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5523 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 5524 5525 /* call Matlab function in ctx with arguments u and y */ 5526 5527 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 5528 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 5529 ierr = PetscMemcpy(&lxdot,&udot,sizeof(u));CHKERRQ(ierr); 5530 ierr = PetscMemcpy(&lA,A,sizeof(u));CHKERRQ(ierr); 5531 ierr = PetscMemcpy(&lB,B,sizeof(u));CHKERRQ(ierr); 5532 5533 prhs[0] = mxCreateDoubleScalar((double)ls); 5534 prhs[1] = mxCreateDoubleScalar((double)time); 5535 prhs[2] = mxCreateDoubleScalar((double)lx); 5536 prhs[3] = mxCreateDoubleScalar((double)lxdot); 5537 prhs[4] = mxCreateDoubleScalar((double)shift); 5538 prhs[5] = mxCreateDoubleScalar((double)lA); 5539 prhs[6] = mxCreateDoubleScalar((double)lB); 5540 prhs[7] = mxCreateString(sctx->funcname); 5541 prhs[8] = sctx->ctx; 5542 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");CHKERRQ(ierr); 5543 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5544 mxDestroyArray(prhs[0]); 5545 mxDestroyArray(prhs[1]); 5546 mxDestroyArray(prhs[2]); 5547 mxDestroyArray(prhs[3]); 5548 mxDestroyArray(prhs[4]); 5549 mxDestroyArray(prhs[5]); 5550 mxDestroyArray(prhs[6]); 5551 mxDestroyArray(prhs[7]); 5552 mxDestroyArray(plhs[0]); 5553 mxDestroyArray(plhs[1]); 5554 PetscFunctionReturn(0); 5555 } 5556 5557 5558 #undef __FUNCT__ 5559 #define __FUNCT__ "TSSetJacobianMatlab" 5560 /* 5561 TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 5562 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 5563 5564 Logically Collective on TS 5565 5566 Input Parameters: 5567 + ts - the TS context 5568 . A,B - Jacobian matrices 5569 . func - function evaluation routine 5570 - ctx - user context 5571 5572 Calling sequence of func: 5573 $ flag = func (TS ts,PetscReal time,Vec u,Vec udot,Mat A,Mat B,void *ctx); 5574 5575 5576 Level: developer 5577 5578 .keywords: TS, nonlinear, set, function 5579 5580 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 5581 */ 5582 PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx) 5583 { 5584 PetscErrorCode ierr; 5585 TSMatlabContext *sctx; 5586 5587 PetscFunctionBegin; 5588 /* currently sctx is memory bleed */ 5589 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 5590 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5591 /* 5592 This should work, but it doesn't 5593 sctx->ctx = ctx; 5594 mexMakeArrayPersistent(sctx->ctx); 5595 */ 5596 sctx->ctx = mxDuplicateArray(ctx); 5597 5598 ierr = TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 5599 PetscFunctionReturn(0); 5600 } 5601 5602 #undef __FUNCT__ 5603 #define __FUNCT__ "TSMonitor_Matlab" 5604 /* 5605 TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab(). 5606 5607 Collective on TS 5608 5609 .seealso: TSSetFunction(), TSGetFunction() 5610 @*/ 5611 PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec u, void *ctx) 5612 { 5613 PetscErrorCode ierr; 5614 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 5615 int nlhs = 1,nrhs = 6; 5616 mxArray *plhs[1],*prhs[6]; 5617 long long int lx = 0,ls = 0; 5618 5619 PetscFunctionBegin; 5620 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5621 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 5622 5623 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 5624 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 5625 5626 prhs[0] = mxCreateDoubleScalar((double)ls); 5627 prhs[1] = mxCreateDoubleScalar((double)it); 5628 prhs[2] = mxCreateDoubleScalar((double)time); 5629 prhs[3] = mxCreateDoubleScalar((double)lx); 5630 prhs[4] = mxCreateString(sctx->funcname); 5631 prhs[5] = sctx->ctx; 5632 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");CHKERRQ(ierr); 5633 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5634 mxDestroyArray(prhs[0]); 5635 mxDestroyArray(prhs[1]); 5636 mxDestroyArray(prhs[2]); 5637 mxDestroyArray(prhs[3]); 5638 mxDestroyArray(prhs[4]); 5639 mxDestroyArray(plhs[0]); 5640 PetscFunctionReturn(0); 5641 } 5642 5643 5644 #undef __FUNCT__ 5645 #define __FUNCT__ "TSMonitorSetMatlab" 5646 /* 5647 TSMonitorSetMatlab - Sets the monitor function from Matlab 5648 5649 Level: developer 5650 5651 .keywords: TS, nonlinear, set, function 5652 5653 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 5654 */ 5655 PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx) 5656 { 5657 PetscErrorCode ierr; 5658 TSMatlabContext *sctx; 5659 5660 PetscFunctionBegin; 5661 /* currently sctx is memory bleed */ 5662 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 5663 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5664 /* 5665 This should work, but it doesn't 5666 sctx->ctx = ctx; 5667 mexMakeArrayPersistent(sctx->ctx); 5668 */ 5669 sctx->ctx = mxDuplicateArray(ctx); 5670 5671 ierr = TSMonitorSet(ts,TSMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 5672 PetscFunctionReturn(0); 5673 } 5674 #endif 5675 5676 #undef __FUNCT__ 5677 #define __FUNCT__ "TSMonitorLGSolution" 5678 /*@C 5679 TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector 5680 in a time based line graph 5681 5682 Collective on TS 5683 5684 Input Parameters: 5685 + ts - the TS context 5686 . step - current time-step 5687 . ptime - current time 5688 - lg - a line graph object 5689 5690 Options Database: 5691 . -ts_monitor_lg_solution_variables 5692 5693 Level: intermediate 5694 5695 Notes: each process in a parallel run displays its component solutions in a separate window 5696 5697 .keywords: TS, vector, monitor, view 5698 5699 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 5700 @*/ 5701 PetscErrorCode TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 5702 { 5703 PetscErrorCode ierr; 5704 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 5705 const PetscScalar *yy; 5706 PetscInt dim; 5707 Vec v; 5708 5709 PetscFunctionBegin; 5710 if (!step) { 5711 PetscDrawAxis axis; 5712 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 5713 ierr = PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");CHKERRQ(ierr); 5714 if (ctx->names && !ctx->displaynames) { 5715 char **displaynames; 5716 PetscBool flg; 5717 5718 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 5719 ierr = PetscMalloc((dim+1)*sizeof(char*),&displaynames);CHKERRQ(ierr); 5720 ierr = PetscMemzero(displaynames,(dim+1)*sizeof(char*));CHKERRQ(ierr); 5721 ierr = PetscOptionsGetStringArray(((PetscObject)ts)->prefix,"-ts_monitor_lg_solution_variables",displaynames,&dim,&flg);CHKERRQ(ierr); 5722 if (flg) { 5723 ierr = TSMonitorLGCtxSetDisplayVariables(ctx,(const char *const *)displaynames);CHKERRQ(ierr); 5724 } 5725 ierr = PetscStrArrayDestroy(&displaynames);CHKERRQ(ierr); 5726 } 5727 if (ctx->displaynames) { 5728 ierr = PetscDrawLGSetDimension(ctx->lg,ctx->ndisplayvariables);CHKERRQ(ierr); 5729 ierr = PetscDrawLGSetLegend(ctx->lg,(const char *const *)ctx->displaynames);CHKERRQ(ierr); 5730 } else if (ctx->names) { 5731 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 5732 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 5733 ierr = PetscDrawLGSetLegend(ctx->lg,(const char *const *)ctx->names);CHKERRQ(ierr); 5734 } 5735 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 5736 } 5737 if (ctx->transform) { 5738 ierr = (*ctx->transform)(ctx->transformctx,u,&v);CHKERRQ(ierr); 5739 } else { 5740 v = u; 5741 } 5742 ierr = VecGetArrayRead(v,&yy);CHKERRQ(ierr); 5743 #if defined(PETSC_USE_COMPLEX) 5744 { 5745 PetscReal *yreal; 5746 PetscInt i,n; 5747 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 5748 ierr = PetscMalloc1(n,&yreal);CHKERRQ(ierr); 5749 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 5750 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 5751 ierr = PetscFree(yreal);CHKERRQ(ierr); 5752 } 5753 #else 5754 if (ctx->displaynames) { 5755 PetscInt i; 5756 for (i=0; i<ctx->ndisplayvariables; i++) { 5757 ctx->displayvalues[i] = yy[ctx->displayvariables[i]]; 5758 } 5759 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,ctx->displayvalues);CHKERRQ(ierr); 5760 } else { 5761 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 5762 } 5763 #endif 5764 ierr = VecRestoreArrayRead(v,&yy);CHKERRQ(ierr); 5765 if (ctx->transform) { 5766 ierr = VecDestroy(&v);CHKERRQ(ierr); 5767 } 5768 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 5769 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 5770 } 5771 PetscFunctionReturn(0); 5772 } 5773 5774 5775 #undef __FUNCT__ 5776 #define __FUNCT__ "TSMonitorLGSetVariableNames" 5777 /*@C 5778 TSMonitorLGSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot 5779 5780 Collective on TS 5781 5782 Input Parameters: 5783 + ts - the TS context 5784 - names - the names of the components, final string must be NULL 5785 5786 Level: intermediate 5787 5788 .keywords: TS, vector, monitor, view 5789 5790 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetVariableNames() 5791 @*/ 5792 PetscErrorCode TSMonitorLGSetVariableNames(TS ts,const char * const *names) 5793 { 5794 PetscErrorCode ierr; 5795 PetscInt i; 5796 5797 PetscFunctionBegin; 5798 for (i=0; i<ts->numbermonitors; i++) { 5799 if (ts->monitor[i] == TSMonitorLGSolution) { 5800 ierr = TSMonitorLGCtxSetVariableNames((TSMonitorLGCtx)ts->monitorcontext[i],names);CHKERRQ(ierr); 5801 break; 5802 } 5803 } 5804 PetscFunctionReturn(0); 5805 } 5806 5807 #undef __FUNCT__ 5808 #define __FUNCT__ "TSMonitorLGCtxSetVariableNames" 5809 /*@C 5810 TSMonitorLGCtxSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot 5811 5812 Collective on TS 5813 5814 Input Parameters: 5815 + ts - the TS context 5816 - names - the names of the components, final string must be NULL 5817 5818 Level: intermediate 5819 5820 .keywords: TS, vector, monitor, view 5821 5822 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables(), TSMonitorLGSetVariableNames() 5823 @*/ 5824 PetscErrorCode TSMonitorLGCtxSetVariableNames(TSMonitorLGCtx ctx,const char * const *names) 5825 { 5826 PetscErrorCode ierr; 5827 5828 PetscFunctionBegin; 5829 ierr = PetscStrArrayDestroy(&ctx->names);CHKERRQ(ierr); 5830 ierr = PetscStrArrayallocpy(names,&ctx->names);CHKERRQ(ierr); 5831 PetscFunctionReturn(0); 5832 } 5833 5834 #undef __FUNCT__ 5835 #define __FUNCT__ "TSMonitorLGGetVariableNames" 5836 /*@C 5837 TSMonitorLGGetVariableNames - Gets the name of each component in the solution vector so that it may be displayed in the plot 5838 5839 Collective on TS 5840 5841 Input Parameter: 5842 . ts - the TS context 5843 5844 Output Parameter: 5845 . names - the names of the components, final string must be NULL 5846 5847 Level: intermediate 5848 5849 .keywords: TS, vector, monitor, view 5850 5851 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables() 5852 @*/ 5853 PetscErrorCode TSMonitorLGGetVariableNames(TS ts,const char *const **names) 5854 { 5855 PetscInt i; 5856 5857 PetscFunctionBegin; 5858 *names = NULL; 5859 for (i=0; i<ts->numbermonitors; i++) { 5860 if (ts->monitor[i] == TSMonitorLGSolution) { 5861 TSMonitorLGCtx ctx = (TSMonitorLGCtx) ts->monitorcontext[i]; 5862 *names = (const char *const *)ctx->names; 5863 break; 5864 } 5865 } 5866 PetscFunctionReturn(0); 5867 } 5868 5869 #undef __FUNCT__ 5870 #define __FUNCT__ "TSMonitorLGCtxSetDisplayVariables" 5871 /*@C 5872 TSMonitorLGCtxSetDisplayVariables - Sets the variables that are to be display in the monitor 5873 5874 Collective on TS 5875 5876 Input Parameters: 5877 + ctx - the TSMonitorLG context 5878 . displaynames - the names of the components, final string must be NULL 5879 5880 Level: intermediate 5881 5882 .keywords: TS, vector, monitor, view 5883 5884 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames() 5885 @*/ 5886 PetscErrorCode TSMonitorLGCtxSetDisplayVariables(TSMonitorLGCtx ctx,const char * const *displaynames) 5887 { 5888 PetscInt j = 0,k; 5889 PetscErrorCode ierr; 5890 5891 PetscFunctionBegin; 5892 if (!ctx->names) PetscFunctionReturn(0); 5893 ierr = PetscStrArrayDestroy(&ctx->displaynames);CHKERRQ(ierr); 5894 ierr = PetscStrArrayallocpy(displaynames,&ctx->displaynames);CHKERRQ(ierr); 5895 while (displaynames[j]) j++; 5896 ctx->ndisplayvariables = j; 5897 ierr = PetscMalloc1(ctx->ndisplayvariables,&ctx->displayvariables);CHKERRQ(ierr); 5898 ierr = PetscMalloc1(ctx->ndisplayvariables,&ctx->displayvalues);CHKERRQ(ierr); 5899 j = 0; 5900 while (displaynames[j]) { 5901 k = 0; 5902 while (ctx->names[k]) { 5903 PetscBool flg; 5904 ierr = PetscStrcmp(displaynames[j],ctx->names[k],&flg);CHKERRQ(ierr); 5905 if (flg) { 5906 ctx->displayvariables[j] = k; 5907 break; 5908 } 5909 k++; 5910 } 5911 j++; 5912 } 5913 PetscFunctionReturn(0); 5914 } 5915 5916 5917 #undef __FUNCT__ 5918 #define __FUNCT__ "TSMonitorLGSetDisplayVariables" 5919 /*@C 5920 TSMonitorLGSetDisplayVariables - Sets the variables that are to be display in the monitor 5921 5922 Collective on TS 5923 5924 Input Parameters: 5925 + ts - the TS context 5926 . displaynames - the names of the components, final string must be NULL 5927 5928 Level: intermediate 5929 5930 .keywords: TS, vector, monitor, view 5931 5932 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames() 5933 @*/ 5934 PetscErrorCode TSMonitorLGSetDisplayVariables(TS ts,const char * const *displaynames) 5935 { 5936 PetscInt i; 5937 PetscErrorCode ierr; 5938 5939 PetscFunctionBegin; 5940 for (i=0; i<ts->numbermonitors; i++) { 5941 if (ts->monitor[i] == TSMonitorLGSolution) { 5942 ierr = TSMonitorLGCtxSetDisplayVariables((TSMonitorLGCtx)ts->monitorcontext[i],displaynames);CHKERRQ(ierr); 5943 break; 5944 } 5945 } 5946 PetscFunctionReturn(0); 5947 } 5948 5949 #undef __FUNCT__ 5950 #define __FUNCT__ "TSMonitorLGSetTransform" 5951 /*@C 5952 TSMonitorLGSetTransform - Solution vector will be transformed by provided function before being displayed 5953 5954 Collective on TS 5955 5956 Input Parameters: 5957 + ts - the TS context 5958 . transform - the transform function 5959 . destroy - function to destroy the optional context 5960 - ctx - optional context used by transform function 5961 5962 Level: intermediate 5963 5964 .keywords: TS, vector, monitor, view 5965 5966 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames(), TSMonitorLGCtxSetTransform() 5967 @*/ 5968 PetscErrorCode TSMonitorLGSetTransform(TS ts,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx) 5969 { 5970 PetscInt i; 5971 PetscErrorCode ierr; 5972 5973 PetscFunctionBegin; 5974 for (i=0; i<ts->numbermonitors; i++) { 5975 if (ts->monitor[i] == TSMonitorLGSolution) { 5976 ierr = TSMonitorLGCtxSetTransform((TSMonitorLGCtx)ts->monitorcontext[i],transform,destroy,tctx);CHKERRQ(ierr); 5977 } 5978 } 5979 PetscFunctionReturn(0); 5980 } 5981 5982 #undef __FUNCT__ 5983 #define __FUNCT__ "TSMonitorLGCtxSetTransform" 5984 /*@C 5985 TSMonitorLGCtxSetTransform - Solution vector will be transformed by provided function before being displayed 5986 5987 Collective on TSLGCtx 5988 5989 Input Parameters: 5990 + ts - the TS context 5991 . transform - the transform function 5992 . destroy - function to destroy the optional context 5993 - ctx - optional context used by transform function 5994 5995 Level: intermediate 5996 5997 .keywords: TS, vector, monitor, view 5998 5999 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames(), TSMonitorLGSetTransform() 6000 @*/ 6001 PetscErrorCode TSMonitorLGCtxSetTransform(TSMonitorLGCtx ctx,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx) 6002 { 6003 PetscFunctionBegin; 6004 ctx->transform = transform; 6005 ctx->transformdestroy = destroy; 6006 ctx->transformctx = tctx; 6007 PetscFunctionReturn(0); 6008 } 6009 6010 #undef __FUNCT__ 6011 #define __FUNCT__ "TSMonitorLGError" 6012 /*@C 6013 TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the solution vector 6014 in a time based line graph 6015 6016 Collective on TS 6017 6018 Input Parameters: 6019 + ts - the TS context 6020 . step - current time-step 6021 . ptime - current time 6022 - lg - a line graph object 6023 6024 Level: intermediate 6025 6026 Notes: 6027 Only for sequential solves. 6028 6029 The user must provide the solution using TSSetSolutionFunction() to use this monitor. 6030 6031 Options Database Keys: 6032 . -ts_monitor_lg_error - create a graphical monitor of error history 6033 6034 .keywords: TS, vector, monitor, view 6035 6036 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction() 6037 @*/ 6038 PetscErrorCode TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 6039 { 6040 PetscErrorCode ierr; 6041 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 6042 const PetscScalar *yy; 6043 Vec y; 6044 PetscInt dim; 6045 6046 PetscFunctionBegin; 6047 if (!step) { 6048 PetscDrawAxis axis; 6049 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 6050 ierr = PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Solution");CHKERRQ(ierr); 6051 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 6052 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 6053 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 6054 } 6055 ierr = VecDuplicate(u,&y);CHKERRQ(ierr); 6056 ierr = TSComputeSolutionFunction(ts,ptime,y);CHKERRQ(ierr); 6057 ierr = VecAXPY(y,-1.0,u);CHKERRQ(ierr); 6058 ierr = VecGetArrayRead(y,&yy);CHKERRQ(ierr); 6059 #if defined(PETSC_USE_COMPLEX) 6060 { 6061 PetscReal *yreal; 6062 PetscInt i,n; 6063 ierr = VecGetLocalSize(y,&n);CHKERRQ(ierr); 6064 ierr = PetscMalloc1(n,&yreal);CHKERRQ(ierr); 6065 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 6066 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 6067 ierr = PetscFree(yreal);CHKERRQ(ierr); 6068 } 6069 #else 6070 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 6071 #endif 6072 ierr = VecRestoreArrayRead(y,&yy);CHKERRQ(ierr); 6073 ierr = VecDestroy(&y);CHKERRQ(ierr); 6074 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 6075 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 6076 } 6077 PetscFunctionReturn(0); 6078 } 6079 6080 #undef __FUNCT__ 6081 #define __FUNCT__ "TSMonitorLGSNESIterations" 6082 PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 6083 { 6084 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 6085 PetscReal x = ptime,y; 6086 PetscErrorCode ierr; 6087 PetscInt its; 6088 6089 PetscFunctionBegin; 6090 if (!n) { 6091 PetscDrawAxis axis; 6092 6093 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 6094 ierr = PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");CHKERRQ(ierr); 6095 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 6096 6097 ctx->snes_its = 0; 6098 } 6099 ierr = TSGetSNESIterations(ts,&its);CHKERRQ(ierr); 6100 y = its - ctx->snes_its; 6101 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 6102 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 6103 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 6104 } 6105 ctx->snes_its = its; 6106 PetscFunctionReturn(0); 6107 } 6108 6109 #undef __FUNCT__ 6110 #define __FUNCT__ "TSMonitorLGKSPIterations" 6111 PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 6112 { 6113 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 6114 PetscReal x = ptime,y; 6115 PetscErrorCode ierr; 6116 PetscInt its; 6117 6118 PetscFunctionBegin; 6119 if (!n) { 6120 PetscDrawAxis axis; 6121 6122 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 6123 ierr = PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");CHKERRQ(ierr); 6124 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 6125 6126 ctx->ksp_its = 0; 6127 } 6128 ierr = TSGetKSPIterations(ts,&its);CHKERRQ(ierr); 6129 y = its - ctx->ksp_its; 6130 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 6131 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 6132 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 6133 } 6134 ctx->ksp_its = its; 6135 PetscFunctionReturn(0); 6136 } 6137 6138 #undef __FUNCT__ 6139 #define __FUNCT__ "TSComputeLinearStability" 6140 /*@ 6141 TSComputeLinearStability - computes the linear stability function at a point 6142 6143 Collective on TS and Vec 6144 6145 Input Parameters: 6146 + ts - the TS context 6147 - xr,xi - real and imaginary part of input arguments 6148 6149 Output Parameters: 6150 . yr,yi - real and imaginary part of function value 6151 6152 Level: developer 6153 6154 .keywords: TS, compute 6155 6156 .seealso: TSSetRHSFunction(), TSComputeIFunction() 6157 @*/ 6158 PetscErrorCode TSComputeLinearStability(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi) 6159 { 6160 PetscErrorCode ierr; 6161 6162 PetscFunctionBegin; 6163 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 6164 if (!ts->ops->linearstability) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Linearized stability function not provided for this method"); 6165 ierr = (*ts->ops->linearstability)(ts,xr,xi,yr,yi);CHKERRQ(ierr); 6166 PetscFunctionReturn(0); 6167 } 6168 6169 /* ------------------------------------------------------------------------*/ 6170 #undef __FUNCT__ 6171 #define __FUNCT__ "TSMonitorEnvelopeCtxCreate" 6172 /*@C 6173 TSMonitorEnvelopeCtxCreate - Creates a context for use with TSMonitorEnvelope() 6174 6175 Collective on TS 6176 6177 Input Parameters: 6178 . ts - the ODE solver object 6179 6180 Output Parameter: 6181 . ctx - the context 6182 6183 Level: intermediate 6184 6185 .keywords: TS, monitor, line graph, residual, seealso 6186 6187 .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError() 6188 6189 @*/ 6190 PetscErrorCode TSMonitorEnvelopeCtxCreate(TS ts,TSMonitorEnvelopeCtx *ctx) 6191 { 6192 PetscErrorCode ierr; 6193 6194 PetscFunctionBegin; 6195 ierr = PetscNew(ctx);CHKERRQ(ierr); 6196 PetscFunctionReturn(0); 6197 } 6198 6199 #undef __FUNCT__ 6200 #define __FUNCT__ "TSMonitorEnvelope" 6201 /*@C 6202 TSMonitorEnvelope - Monitors the maximum and minimum value of each component of the solution 6203 6204 Collective on TS 6205 6206 Input Parameters: 6207 + ts - the TS context 6208 . step - current time-step 6209 . ptime - current time 6210 - ctx - the envelope context 6211 6212 Options Database: 6213 . -ts_monitor_envelope 6214 6215 Level: intermediate 6216 6217 Notes: after a solve you can use TSMonitorEnvelopeGetBounds() to access the envelope 6218 6219 .keywords: TS, vector, monitor, view 6220 6221 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorEnvelopeGetBounds() 6222 @*/ 6223 PetscErrorCode TSMonitorEnvelope(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 6224 { 6225 PetscErrorCode ierr; 6226 TSMonitorEnvelopeCtx ctx = (TSMonitorEnvelopeCtx)dummy; 6227 6228 PetscFunctionBegin; 6229 if (!ctx->max) { 6230 ierr = VecDuplicate(u,&ctx->max);CHKERRQ(ierr); 6231 ierr = VecDuplicate(u,&ctx->min);CHKERRQ(ierr); 6232 ierr = VecCopy(u,ctx->max);CHKERRQ(ierr); 6233 ierr = VecCopy(u,ctx->min);CHKERRQ(ierr); 6234 } else { 6235 ierr = VecPointwiseMax(ctx->max,u,ctx->max);CHKERRQ(ierr); 6236 ierr = VecPointwiseMin(ctx->min,u,ctx->min);CHKERRQ(ierr); 6237 } 6238 PetscFunctionReturn(0); 6239 } 6240 6241 6242 #undef __FUNCT__ 6243 #define __FUNCT__ "TSMonitorEnvelopeGetBounds" 6244 /*@C 6245 TSMonitorEnvelopeGetBounds - Gets the bounds for the components of the solution 6246 6247 Collective on TS 6248 6249 Input Parameter: 6250 . ts - the TS context 6251 6252 Output Parameter: 6253 + max - the maximum values 6254 - min - the minimum values 6255 6256 Level: intermediate 6257 6258 .keywords: TS, vector, monitor, view 6259 6260 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables() 6261 @*/ 6262 PetscErrorCode TSMonitorEnvelopeGetBounds(TS ts,Vec *max,Vec *min) 6263 { 6264 PetscInt i; 6265 6266 PetscFunctionBegin; 6267 if (max) *max = NULL; 6268 if (min) *min = NULL; 6269 for (i=0; i<ts->numbermonitors; i++) { 6270 if (ts->monitor[i] == TSMonitorEnvelope) { 6271 TSMonitorEnvelopeCtx ctx = (TSMonitorEnvelopeCtx) ts->monitorcontext[i]; 6272 if (max) *max = ctx->max; 6273 if (min) *min = ctx->min; 6274 break; 6275 } 6276 } 6277 PetscFunctionReturn(0); 6278 } 6279 6280 #undef __FUNCT__ 6281 #define __FUNCT__ "TSMonitorEnvelopeCtxDestroy" 6282 /*@C 6283 TSMonitorEnvelopeCtxDestroy - Destroys a context that was created with TSMonitorEnvelopeCtxCreate(). 6284 6285 Collective on TSMonitorEnvelopeCtx 6286 6287 Input Parameter: 6288 . ctx - the monitor context 6289 6290 Level: intermediate 6291 6292 .keywords: TS, monitor, line graph, destroy 6293 6294 .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep(); 6295 @*/ 6296 PetscErrorCode TSMonitorEnvelopeCtxDestroy(TSMonitorEnvelopeCtx *ctx) 6297 { 6298 PetscErrorCode ierr; 6299 6300 PetscFunctionBegin; 6301 ierr = VecDestroy(&(*ctx)->min);CHKERRQ(ierr); 6302 ierr = VecDestroy(&(*ctx)->max);CHKERRQ(ierr); 6303 ierr = PetscFree(*ctx);CHKERRQ(ierr); 6304 PetscFunctionReturn(0); 6305 } 6306 6307 #undef __FUNCT__ 6308 #define __FUNCT__ "TSRollBack" 6309 /*@ 6310 TSRollBack - Rolls back one time step 6311 6312 Collective on TS 6313 6314 Input Parameter: 6315 . ts - the TS context obtained from TSCreate() 6316 6317 Level: advanced 6318 6319 .keywords: TS, timestep, rollback 6320 6321 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSInterpolate() 6322 @*/ 6323 PetscErrorCode TSRollBack(TS ts) 6324 { 6325 PetscErrorCode ierr; 6326 6327 PetscFunctionBegin; 6328 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 6329 6330 if (!ts->ops->rollback) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRollBack not implemented for type '%s'",((PetscObject)ts)->type_name); 6331 ierr = (*ts->ops->rollback)(ts);CHKERRQ(ierr); 6332 ts->time_step = ts->ptime - ts->ptime_prev; 6333 ts->ptime = ts->ptime_prev; 6334 ts->steprollback = PETSC_TRUE; /* Flag to indicate that the step is rollbacked */ 6335 PetscFunctionReturn(0); 6336 } 6337 6338 #undef __FUNCT__ 6339 #define __FUNCT__ "TSGetStages" 6340 /*@ 6341 TSGetStages - Get the number of stages and stage values 6342 6343 Input Parameter: 6344 . ts - the TS context obtained from TSCreate() 6345 6346 Level: advanced 6347 6348 .keywords: TS, getstages 6349 6350 .seealso: TSCreate() 6351 @*/ 6352 PetscErrorCode TSGetStages(TS ts,PetscInt *ns, Vec **Y) 6353 { 6354 PetscErrorCode ierr; 6355 6356 PetscFunctionBegin; 6357 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 6358 PetscValidPointer(ns,2); 6359 6360 if (!ts->ops->getstages) *ns=0; 6361 else { 6362 ierr = (*ts->ops->getstages)(ts,ns,Y);CHKERRQ(ierr); 6363 } 6364 PetscFunctionReturn(0); 6365 } 6366 6367