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