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 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 PetscFunctionReturn(0); 4051 } 4052 4053 #undef __FUNCT__ 4054 #define __FUNCT__ "TSMonitorDrawSolutionPhase" 4055 /*@C 4056 TSMonitorDrawSolutionPhase - Monitors progress of the TS solvers by plotting the solution as a phase diagram 4057 4058 Collective on TS 4059 4060 Input Parameters: 4061 + ts - the TS context 4062 . step - current time-step 4063 . ptime - current time 4064 - dummy - either a viewer or NULL 4065 4066 Level: intermediate 4067 4068 .keywords: TS, vector, monitor, view 4069 4070 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4071 @*/ 4072 PetscErrorCode TSMonitorDrawSolutionPhase(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4073 { 4074 PetscErrorCode ierr; 4075 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy; 4076 PetscDraw draw; 4077 MPI_Comm comm; 4078 PetscInt n; 4079 PetscMPIInt size; 4080 PetscReal xl,yl,xr,yr,h; 4081 char time[32]; 4082 const PetscScalar *U; 4083 4084 PetscFunctionBegin; 4085 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 4086 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4087 if (size != 1) SETERRQ(comm,PETSC_ERR_SUP,"Only allowed for sequential runs"); 4088 ierr = VecGetSize(u,&n);CHKERRQ(ierr); 4089 if (n != 2) SETERRQ(comm,PETSC_ERR_SUP,"Only for ODEs with two unknowns"); 4090 4091 ierr = PetscViewerDrawGetDraw(ictx->viewer,0,&draw);CHKERRQ(ierr); 4092 4093 ierr = VecGetArrayRead(u,&U);CHKERRQ(ierr); 4094 ierr = PetscDrawAxisGetLimits(ictx->axis,&xl,&xr,&yl,&yr);CHKERRQ(ierr); 4095 if ((PetscRealPart(U[0]) < xl) || (PetscRealPart(U[1]) < yl) || (PetscRealPart(U[0]) > xr) || (PetscRealPart(U[1]) > yr)) { 4096 ierr = VecRestoreArrayRead(u,&U);CHKERRQ(ierr); 4097 PetscFunctionReturn(0); 4098 } 4099 if (!step) ictx->color++; 4100 ierr = PetscDrawPoint(draw,PetscRealPart(U[0]),PetscRealPart(U[1]),ictx->color);CHKERRQ(ierr); 4101 ierr = VecRestoreArrayRead(u,&U);CHKERRQ(ierr); 4102 4103 if (ictx->showtimestepandtime) { 4104 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 4105 ierr = PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);CHKERRQ(ierr); 4106 h = yl + .95*(yr - yl); 4107 ierr = PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);CHKERRQ(ierr); 4108 } 4109 ierr = PetscDrawFlush(draw);CHKERRQ(ierr); 4110 PetscFunctionReturn(0); 4111 } 4112 4113 4114 #undef __FUNCT__ 4115 #define __FUNCT__ "TSMonitorDrawCtxDestroy" 4116 /*@C 4117 TSMonitorDrawCtxDestroy - Destroys the monitor context for TSMonitorDrawSolution() 4118 4119 Collective on TS 4120 4121 Input Parameters: 4122 . ctx - the monitor context 4123 4124 Level: intermediate 4125 4126 .keywords: TS, vector, monitor, view 4127 4128 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawSolution(), TSMonitorDrawError() 4129 @*/ 4130 PetscErrorCode TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx) 4131 { 4132 PetscErrorCode ierr; 4133 4134 PetscFunctionBegin; 4135 ierr = PetscDrawAxisDestroy(&(*ictx)->axis);CHKERRQ(ierr); 4136 ierr = PetscViewerDestroy(&(*ictx)->viewer);CHKERRQ(ierr); 4137 ierr = VecDestroy(&(*ictx)->initialsolution);CHKERRQ(ierr); 4138 ierr = PetscFree(*ictx);CHKERRQ(ierr); 4139 PetscFunctionReturn(0); 4140 } 4141 4142 #undef __FUNCT__ 4143 #define __FUNCT__ "TSMonitorDrawCtxCreate" 4144 /*@C 4145 TSMonitorDrawCtxCreate - Creates the monitor context for TSMonitorDrawCtx 4146 4147 Collective on TS 4148 4149 Input Parameter: 4150 . ts - time-step context 4151 4152 Output Patameter: 4153 . ctx - the monitor context 4154 4155 Options Database: 4156 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 4157 4158 Level: intermediate 4159 4160 .keywords: TS, vector, monitor, view 4161 4162 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawCtx() 4163 @*/ 4164 PetscErrorCode TSMonitorDrawCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorDrawCtx *ctx) 4165 { 4166 PetscErrorCode ierr; 4167 4168 PetscFunctionBegin; 4169 ierr = PetscNew(ctx);CHKERRQ(ierr); 4170 ierr = PetscViewerDrawOpen(comm,host,label,x,y,m,n,&(*ctx)->viewer);CHKERRQ(ierr); 4171 ierr = PetscViewerSetFromOptions((*ctx)->viewer);CHKERRQ(ierr); 4172 4173 (*ctx)->howoften = howoften; 4174 (*ctx)->showinitial = PETSC_FALSE; 4175 ierr = PetscOptionsGetBool(NULL,"-ts_monitor_draw_solution_initial",&(*ctx)->showinitial,NULL);CHKERRQ(ierr); 4176 4177 (*ctx)->showtimestepandtime = PETSC_FALSE; 4178 ierr = PetscOptionsGetBool(NULL,"-ts_monitor_draw_solution_show_time",&(*ctx)->showtimestepandtime,NULL);CHKERRQ(ierr); 4179 (*ctx)->color = PETSC_DRAW_WHITE; 4180 PetscFunctionReturn(0); 4181 } 4182 4183 #undef __FUNCT__ 4184 #define __FUNCT__ "TSMonitorDrawError" 4185 /*@C 4186 TSMonitorDrawError - Monitors progress of the TS solvers by calling 4187 VecView() for the error at each timestep 4188 4189 Collective on TS 4190 4191 Input Parameters: 4192 + ts - the TS context 4193 . step - current time-step 4194 . ptime - current time 4195 - dummy - either a viewer or NULL 4196 4197 Level: intermediate 4198 4199 .keywords: TS, vector, monitor, view 4200 4201 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4202 @*/ 4203 PetscErrorCode TSMonitorDrawError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4204 { 4205 PetscErrorCode ierr; 4206 TSMonitorDrawCtx ctx = (TSMonitorDrawCtx)dummy; 4207 PetscViewer viewer = ctx->viewer; 4208 Vec work; 4209 4210 PetscFunctionBegin; 4211 if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) PetscFunctionReturn(0); 4212 ierr = VecDuplicate(u,&work);CHKERRQ(ierr); 4213 ierr = TSComputeSolutionFunction(ts,ptime,work);CHKERRQ(ierr); 4214 ierr = VecAXPY(work,-1.0,u);CHKERRQ(ierr); 4215 ierr = VecView(work,viewer);CHKERRQ(ierr); 4216 ierr = VecDestroy(&work);CHKERRQ(ierr); 4217 PetscFunctionReturn(0); 4218 } 4219 4220 #include <petsc/private/dmimpl.h> 4221 #undef __FUNCT__ 4222 #define __FUNCT__ "TSSetDM" 4223 /*@ 4224 TSSetDM - Sets the DM that may be used by some preconditioners 4225 4226 Logically Collective on TS and DM 4227 4228 Input Parameters: 4229 + ts - the preconditioner context 4230 - dm - the dm 4231 4232 Level: intermediate 4233 4234 4235 .seealso: TSGetDM(), SNESSetDM(), SNESGetDM() 4236 @*/ 4237 PetscErrorCode TSSetDM(TS ts,DM dm) 4238 { 4239 PetscErrorCode ierr; 4240 SNES snes; 4241 DMTS tsdm; 4242 4243 PetscFunctionBegin; 4244 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4245 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 4246 if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */ 4247 if (ts->dm->dmts && !dm->dmts) { 4248 ierr = DMCopyDMTS(ts->dm,dm);CHKERRQ(ierr); 4249 ierr = DMGetDMTS(ts->dm,&tsdm);CHKERRQ(ierr); 4250 if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */ 4251 tsdm->originaldm = dm; 4252 } 4253 } 4254 ierr = DMDestroy(&ts->dm);CHKERRQ(ierr); 4255 } 4256 ts->dm = dm; 4257 4258 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 4259 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 4260 PetscFunctionReturn(0); 4261 } 4262 4263 #undef __FUNCT__ 4264 #define __FUNCT__ "TSGetDM" 4265 /*@ 4266 TSGetDM - Gets the DM that may be used by some preconditioners 4267 4268 Not Collective 4269 4270 Input Parameter: 4271 . ts - the preconditioner context 4272 4273 Output Parameter: 4274 . dm - the dm 4275 4276 Level: intermediate 4277 4278 4279 .seealso: TSSetDM(), SNESSetDM(), SNESGetDM() 4280 @*/ 4281 PetscErrorCode TSGetDM(TS ts,DM *dm) 4282 { 4283 PetscErrorCode ierr; 4284 4285 PetscFunctionBegin; 4286 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4287 if (!ts->dm) { 4288 ierr = DMShellCreate(PetscObjectComm((PetscObject)ts),&ts->dm);CHKERRQ(ierr); 4289 if (ts->snes) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 4290 } 4291 *dm = ts->dm; 4292 PetscFunctionReturn(0); 4293 } 4294 4295 #undef __FUNCT__ 4296 #define __FUNCT__ "SNESTSFormFunction" 4297 /*@ 4298 SNESTSFormFunction - Function to evaluate nonlinear residual 4299 4300 Logically Collective on SNES 4301 4302 Input Parameter: 4303 + snes - nonlinear solver 4304 . U - the current state at which to evaluate the residual 4305 - ctx - user context, must be a TS 4306 4307 Output Parameter: 4308 . F - the nonlinear residual 4309 4310 Notes: 4311 This function is not normally called by users and is automatically registered with the SNES used by TS. 4312 It is most frequently passed to MatFDColoringSetFunction(). 4313 4314 Level: advanced 4315 4316 .seealso: SNESSetFunction(), MatFDColoringSetFunction() 4317 @*/ 4318 PetscErrorCode SNESTSFormFunction(SNES snes,Vec U,Vec F,void *ctx) 4319 { 4320 TS ts = (TS)ctx; 4321 PetscErrorCode ierr; 4322 4323 PetscFunctionBegin; 4324 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4325 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 4326 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 4327 PetscValidHeaderSpecific(ts,TS_CLASSID,4); 4328 ierr = (ts->ops->snesfunction)(snes,U,F,ts);CHKERRQ(ierr); 4329 PetscFunctionReturn(0); 4330 } 4331 4332 #undef __FUNCT__ 4333 #define __FUNCT__ "SNESTSFormJacobian" 4334 /*@ 4335 SNESTSFormJacobian - Function to evaluate the Jacobian 4336 4337 Collective on SNES 4338 4339 Input Parameter: 4340 + snes - nonlinear solver 4341 . U - the current state at which to evaluate the residual 4342 - ctx - user context, must be a TS 4343 4344 Output Parameter: 4345 + A - the Jacobian 4346 . B - the preconditioning matrix (may be the same as A) 4347 - flag - indicates any structure change in the matrix 4348 4349 Notes: 4350 This function is not normally called by users and is automatically registered with the SNES used by TS. 4351 4352 Level: developer 4353 4354 .seealso: SNESSetJacobian() 4355 @*/ 4356 PetscErrorCode SNESTSFormJacobian(SNES snes,Vec U,Mat A,Mat B,void *ctx) 4357 { 4358 TS ts = (TS)ctx; 4359 PetscErrorCode ierr; 4360 4361 PetscFunctionBegin; 4362 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4363 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 4364 PetscValidPointer(A,3); 4365 PetscValidHeaderSpecific(A,MAT_CLASSID,3); 4366 PetscValidPointer(B,4); 4367 PetscValidHeaderSpecific(B,MAT_CLASSID,4); 4368 PetscValidHeaderSpecific(ts,TS_CLASSID,6); 4369 ierr = (ts->ops->snesjacobian)(snes,U,A,B,ts);CHKERRQ(ierr); 4370 PetscFunctionReturn(0); 4371 } 4372 4373 #undef __FUNCT__ 4374 #define __FUNCT__ "TSComputeRHSFunctionLinear" 4375 /*@C 4376 TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems only 4377 4378 Collective on TS 4379 4380 Input Arguments: 4381 + ts - time stepping context 4382 . t - time at which to evaluate 4383 . U - state at which to evaluate 4384 - ctx - context 4385 4386 Output Arguments: 4387 . F - right hand side 4388 4389 Level: intermediate 4390 4391 Notes: 4392 This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems. 4393 The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian(). 4394 4395 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant() 4396 @*/ 4397 PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec U,Vec F,void *ctx) 4398 { 4399 PetscErrorCode ierr; 4400 Mat Arhs,Brhs; 4401 4402 PetscFunctionBegin; 4403 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 4404 ierr = TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);CHKERRQ(ierr); 4405 ierr = MatMult(Arhs,U,F);CHKERRQ(ierr); 4406 PetscFunctionReturn(0); 4407 } 4408 4409 #undef __FUNCT__ 4410 #define __FUNCT__ "TSComputeRHSJacobianConstant" 4411 /*@C 4412 TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent. 4413 4414 Collective on TS 4415 4416 Input Arguments: 4417 + ts - time stepping context 4418 . t - time at which to evaluate 4419 . U - state at which to evaluate 4420 - ctx - context 4421 4422 Output Arguments: 4423 + A - pointer to operator 4424 . B - pointer to preconditioning matrix 4425 - flg - matrix structure flag 4426 4427 Level: intermediate 4428 4429 Notes: 4430 This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems. 4431 4432 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear() 4433 @*/ 4434 PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx) 4435 { 4436 PetscFunctionBegin; 4437 PetscFunctionReturn(0); 4438 } 4439 4440 #undef __FUNCT__ 4441 #define __FUNCT__ "TSComputeIFunctionLinear" 4442 /*@C 4443 TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only 4444 4445 Collective on TS 4446 4447 Input Arguments: 4448 + ts - time stepping context 4449 . t - time at which to evaluate 4450 . U - state at which to evaluate 4451 . Udot - time derivative of state vector 4452 - ctx - context 4453 4454 Output Arguments: 4455 . F - left hand side 4456 4457 Level: intermediate 4458 4459 Notes: 4460 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 4461 user is required to write their own TSComputeIFunction. 4462 This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems. 4463 The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian(). 4464 4465 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant() 4466 @*/ 4467 PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx) 4468 { 4469 PetscErrorCode ierr; 4470 Mat A,B; 4471 4472 PetscFunctionBegin; 4473 ierr = TSGetIJacobian(ts,&A,&B,NULL,NULL);CHKERRQ(ierr); 4474 ierr = TSComputeIJacobian(ts,t,U,Udot,1.0,A,B,PETSC_TRUE);CHKERRQ(ierr); 4475 ierr = MatMult(A,Udot,F);CHKERRQ(ierr); 4476 PetscFunctionReturn(0); 4477 } 4478 4479 #undef __FUNCT__ 4480 #define __FUNCT__ "TSComputeIJacobianConstant" 4481 /*@C 4482 TSComputeIJacobianConstant - Reuses a time-independent for a semi-implicit DAE or ODE 4483 4484 Collective on TS 4485 4486 Input Arguments: 4487 + ts - time stepping context 4488 . t - time at which to evaluate 4489 . U - state at which to evaluate 4490 . Udot - time derivative of state vector 4491 . shift - shift to apply 4492 - ctx - context 4493 4494 Output Arguments: 4495 + A - pointer to operator 4496 . B - pointer to preconditioning matrix 4497 - flg - matrix structure flag 4498 4499 Level: advanced 4500 4501 Notes: 4502 This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems. 4503 4504 It is only appropriate for problems of the form 4505 4506 $ M Udot = F(U,t) 4507 4508 where M is constant and F is non-stiff. The user must pass M to TSSetIJacobian(). The current implementation only 4509 works with IMEX time integration methods such as TSROSW and TSARKIMEX, since there is no support for de-constructing 4510 an implicit operator of the form 4511 4512 $ shift*M + J 4513 4514 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 4515 a copy of M or reassemble it when requested. 4516 4517 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear() 4518 @*/ 4519 PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat A,Mat B,void *ctx) 4520 { 4521 PetscErrorCode ierr; 4522 4523 PetscFunctionBegin; 4524 ierr = MatScale(A, shift / ts->ijacobian.shift);CHKERRQ(ierr); 4525 ts->ijacobian.shift = shift; 4526 PetscFunctionReturn(0); 4527 } 4528 4529 #undef __FUNCT__ 4530 #define __FUNCT__ "TSGetEquationType" 4531 /*@ 4532 TSGetEquationType - Gets the type of the equation that TS is solving. 4533 4534 Not Collective 4535 4536 Input Parameter: 4537 . ts - the TS context 4538 4539 Output Parameter: 4540 . equation_type - see TSEquationType 4541 4542 Level: beginner 4543 4544 .keywords: TS, equation type 4545 4546 .seealso: TSSetEquationType(), TSEquationType 4547 @*/ 4548 PetscErrorCode TSGetEquationType(TS ts,TSEquationType *equation_type) 4549 { 4550 PetscFunctionBegin; 4551 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4552 PetscValidPointer(equation_type,2); 4553 *equation_type = ts->equation_type; 4554 PetscFunctionReturn(0); 4555 } 4556 4557 #undef __FUNCT__ 4558 #define __FUNCT__ "TSSetEquationType" 4559 /*@ 4560 TSSetEquationType - Sets the type of the equation that TS is solving. 4561 4562 Not Collective 4563 4564 Input Parameter: 4565 + ts - the TS context 4566 - equation_type - see TSEquationType 4567 4568 Level: advanced 4569 4570 .keywords: TS, equation type 4571 4572 .seealso: TSGetEquationType(), TSEquationType 4573 @*/ 4574 PetscErrorCode TSSetEquationType(TS ts,TSEquationType equation_type) 4575 { 4576 PetscFunctionBegin; 4577 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4578 ts->equation_type = equation_type; 4579 PetscFunctionReturn(0); 4580 } 4581 4582 #undef __FUNCT__ 4583 #define __FUNCT__ "TSGetConvergedReason" 4584 /*@ 4585 TSGetConvergedReason - Gets the reason the TS iteration was stopped. 4586 4587 Not Collective 4588 4589 Input Parameter: 4590 . ts - the TS context 4591 4592 Output Parameter: 4593 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 4594 manual pages for the individual convergence tests for complete lists 4595 4596 Level: beginner 4597 4598 Notes: 4599 Can only be called after the call to TSSolve() is complete. 4600 4601 .keywords: TS, nonlinear, set, convergence, test 4602 4603 .seealso: TSSetConvergenceTest(), TSConvergedReason 4604 @*/ 4605 PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason) 4606 { 4607 PetscFunctionBegin; 4608 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4609 PetscValidPointer(reason,2); 4610 *reason = ts->reason; 4611 PetscFunctionReturn(0); 4612 } 4613 4614 #undef __FUNCT__ 4615 #define __FUNCT__ "TSSetConvergedReason" 4616 /*@ 4617 TSSetConvergedReason - Sets the reason for handling the convergence of TSSolve. 4618 4619 Not Collective 4620 4621 Input Parameter: 4622 + ts - the TS context 4623 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 4624 manual pages for the individual convergence tests for complete lists 4625 4626 Level: advanced 4627 4628 Notes: 4629 Can only be called during TSSolve() is active. 4630 4631 .keywords: TS, nonlinear, set, convergence, test 4632 4633 .seealso: TSConvergedReason 4634 @*/ 4635 PetscErrorCode TSSetConvergedReason(TS ts,TSConvergedReason reason) 4636 { 4637 PetscFunctionBegin; 4638 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4639 ts->reason = reason; 4640 PetscFunctionReturn(0); 4641 } 4642 4643 #undef __FUNCT__ 4644 #define __FUNCT__ "TSGetSolveTime" 4645 /*@ 4646 TSGetSolveTime - Gets the time after a call to TSSolve() 4647 4648 Not Collective 4649 4650 Input Parameter: 4651 . ts - the TS context 4652 4653 Output Parameter: 4654 . ftime - the final time. This time should correspond to the final time set with TSSetDuration() 4655 4656 Level: beginner 4657 4658 Notes: 4659 Can only be called after the call to TSSolve() is complete. 4660 4661 .keywords: TS, nonlinear, set, convergence, test 4662 4663 .seealso: TSSetConvergenceTest(), TSConvergedReason 4664 @*/ 4665 PetscErrorCode TSGetSolveTime(TS ts,PetscReal *ftime) 4666 { 4667 PetscFunctionBegin; 4668 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4669 PetscValidPointer(ftime,2); 4670 *ftime = ts->solvetime; 4671 PetscFunctionReturn(0); 4672 } 4673 4674 #undef __FUNCT__ 4675 #define __FUNCT__ "TSGetTotalSteps" 4676 /*@ 4677 TSGetTotalSteps - Gets the total number of steps done since the last call to TSSetUp() or TSCreate() 4678 4679 Not Collective 4680 4681 Input Parameter: 4682 . ts - the TS context 4683 4684 Output Parameter: 4685 . steps - the number of steps 4686 4687 Level: beginner 4688 4689 Notes: 4690 Includes the number of steps for all calls to TSSolve() since TSSetUp() was called 4691 4692 .keywords: TS, nonlinear, set, convergence, test 4693 4694 .seealso: TSSetConvergenceTest(), TSConvergedReason 4695 @*/ 4696 PetscErrorCode TSGetTotalSteps(TS ts,PetscInt *steps) 4697 { 4698 PetscFunctionBegin; 4699 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4700 PetscValidPointer(steps,2); 4701 *steps = ts->total_steps; 4702 PetscFunctionReturn(0); 4703 } 4704 4705 #undef __FUNCT__ 4706 #define __FUNCT__ "TSGetSNESIterations" 4707 /*@ 4708 TSGetSNESIterations - Gets the total number of nonlinear iterations 4709 used by the time integrator. 4710 4711 Not Collective 4712 4713 Input Parameter: 4714 . ts - TS context 4715 4716 Output Parameter: 4717 . nits - number of nonlinear iterations 4718 4719 Notes: 4720 This counter is reset to zero for each successive call to TSSolve(). 4721 4722 Level: intermediate 4723 4724 .keywords: TS, get, number, nonlinear, iterations 4725 4726 .seealso: TSGetKSPIterations() 4727 @*/ 4728 PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits) 4729 { 4730 PetscFunctionBegin; 4731 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4732 PetscValidIntPointer(nits,2); 4733 *nits = ts->snes_its; 4734 PetscFunctionReturn(0); 4735 } 4736 4737 #undef __FUNCT__ 4738 #define __FUNCT__ "TSGetKSPIterations" 4739 /*@ 4740 TSGetKSPIterations - Gets the total number of linear iterations 4741 used by the time integrator. 4742 4743 Not Collective 4744 4745 Input Parameter: 4746 . ts - TS context 4747 4748 Output Parameter: 4749 . lits - number of linear iterations 4750 4751 Notes: 4752 This counter is reset to zero for each successive call to TSSolve(). 4753 4754 Level: intermediate 4755 4756 .keywords: TS, get, number, linear, iterations 4757 4758 .seealso: TSGetSNESIterations(), SNESGetKSPIterations() 4759 @*/ 4760 PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits) 4761 { 4762 PetscFunctionBegin; 4763 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4764 PetscValidIntPointer(lits,2); 4765 *lits = ts->ksp_its; 4766 PetscFunctionReturn(0); 4767 } 4768 4769 #undef __FUNCT__ 4770 #define __FUNCT__ "TSGetStepRejections" 4771 /*@ 4772 TSGetStepRejections - Gets the total number of rejected steps. 4773 4774 Not Collective 4775 4776 Input Parameter: 4777 . ts - TS context 4778 4779 Output Parameter: 4780 . rejects - number of steps rejected 4781 4782 Notes: 4783 This counter is reset to zero for each successive call to TSSolve(). 4784 4785 Level: intermediate 4786 4787 .keywords: TS, get, number 4788 4789 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails() 4790 @*/ 4791 PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects) 4792 { 4793 PetscFunctionBegin; 4794 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4795 PetscValidIntPointer(rejects,2); 4796 *rejects = ts->reject; 4797 PetscFunctionReturn(0); 4798 } 4799 4800 #undef __FUNCT__ 4801 #define __FUNCT__ "TSGetSNESFailures" 4802 /*@ 4803 TSGetSNESFailures - Gets the total number of failed SNES solves 4804 4805 Not Collective 4806 4807 Input Parameter: 4808 . ts - TS context 4809 4810 Output Parameter: 4811 . fails - number of failed nonlinear solves 4812 4813 Notes: 4814 This counter is reset to zero for each successive call to TSSolve(). 4815 4816 Level: intermediate 4817 4818 .keywords: TS, get, number 4819 4820 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures() 4821 @*/ 4822 PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails) 4823 { 4824 PetscFunctionBegin; 4825 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4826 PetscValidIntPointer(fails,2); 4827 *fails = ts->num_snes_failures; 4828 PetscFunctionReturn(0); 4829 } 4830 4831 #undef __FUNCT__ 4832 #define __FUNCT__ "TSSetMaxStepRejections" 4833 /*@ 4834 TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails 4835 4836 Not Collective 4837 4838 Input Parameter: 4839 + ts - TS context 4840 - rejects - maximum number of rejected steps, pass -1 for unlimited 4841 4842 Notes: 4843 The counter is reset to zero for each step 4844 4845 Options Database Key: 4846 . -ts_max_reject - Maximum number of step rejections before a step fails 4847 4848 Level: intermediate 4849 4850 .keywords: TS, set, maximum, number 4851 4852 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 4853 @*/ 4854 PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects) 4855 { 4856 PetscFunctionBegin; 4857 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4858 ts->max_reject = rejects; 4859 PetscFunctionReturn(0); 4860 } 4861 4862 #undef __FUNCT__ 4863 #define __FUNCT__ "TSSetMaxSNESFailures" 4864 /*@ 4865 TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves 4866 4867 Not Collective 4868 4869 Input Parameter: 4870 + ts - TS context 4871 - fails - maximum number of failed nonlinear solves, pass -1 for unlimited 4872 4873 Notes: 4874 The counter is reset to zero for each successive call to TSSolve(). 4875 4876 Options Database Key: 4877 . -ts_max_snes_failures - Maximum number of nonlinear solve failures 4878 4879 Level: intermediate 4880 4881 .keywords: TS, set, maximum, number 4882 4883 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason() 4884 @*/ 4885 PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails) 4886 { 4887 PetscFunctionBegin; 4888 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4889 ts->max_snes_failures = fails; 4890 PetscFunctionReturn(0); 4891 } 4892 4893 #undef __FUNCT__ 4894 #define __FUNCT__ "TSSetErrorIfStepFails" 4895 /*@ 4896 TSSetErrorIfStepFails - Error if no step succeeds 4897 4898 Not Collective 4899 4900 Input Parameter: 4901 + ts - TS context 4902 - err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure 4903 4904 Options Database Key: 4905 . -ts_error_if_step_fails - Error if no step succeeds 4906 4907 Level: intermediate 4908 4909 .keywords: TS, set, error 4910 4911 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 4912 @*/ 4913 PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err) 4914 { 4915 PetscFunctionBegin; 4916 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4917 ts->errorifstepfailed = err; 4918 PetscFunctionReturn(0); 4919 } 4920 4921 #undef __FUNCT__ 4922 #define __FUNCT__ "TSMonitorSolutionBinary" 4923 /*@C 4924 TSMonitorSolutionBinary - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file 4925 4926 Collective on TS 4927 4928 Input Parameters: 4929 + ts - the TS context 4930 . step - current time-step 4931 . ptime - current time 4932 . u - current state 4933 - viewer - binary viewer 4934 4935 Level: intermediate 4936 4937 .keywords: TS, vector, monitor, view 4938 4939 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4940 @*/ 4941 PetscErrorCode TSMonitorSolutionBinary(TS ts,PetscInt step,PetscReal ptime,Vec u,void *viewer) 4942 { 4943 PetscErrorCode ierr; 4944 PetscViewer v = (PetscViewer)viewer; 4945 4946 PetscFunctionBegin; 4947 ierr = VecView(u,v);CHKERRQ(ierr); 4948 PetscFunctionReturn(0); 4949 } 4950 4951 #undef __FUNCT__ 4952 #define __FUNCT__ "TSMonitorSolutionVTK" 4953 /*@C 4954 TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep. 4955 4956 Collective on TS 4957 4958 Input Parameters: 4959 + ts - the TS context 4960 . step - current time-step 4961 . ptime - current time 4962 . u - current state 4963 - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 4964 4965 Level: intermediate 4966 4967 Notes: 4968 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. 4969 These are named according to the file name template. 4970 4971 This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy(). 4972 4973 .keywords: TS, vector, monitor, view 4974 4975 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4976 @*/ 4977 PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec u,void *filenametemplate) 4978 { 4979 PetscErrorCode ierr; 4980 char filename[PETSC_MAX_PATH_LEN]; 4981 PetscViewer viewer; 4982 4983 PetscFunctionBegin; 4984 ierr = PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);CHKERRQ(ierr); 4985 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)ts),filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 4986 ierr = VecView(u,viewer);CHKERRQ(ierr); 4987 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 4988 PetscFunctionReturn(0); 4989 } 4990 4991 #undef __FUNCT__ 4992 #define __FUNCT__ "TSMonitorSolutionVTKDestroy" 4993 /*@C 4994 TSMonitorSolutionVTKDestroy - Destroy context for monitoring 4995 4996 Collective on TS 4997 4998 Input Parameters: 4999 . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 5000 5001 Level: intermediate 5002 5003 Note: 5004 This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK(). 5005 5006 .keywords: TS, vector, monitor, view 5007 5008 .seealso: TSMonitorSet(), TSMonitorSolutionVTK() 5009 @*/ 5010 PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate) 5011 { 5012 PetscErrorCode ierr; 5013 5014 PetscFunctionBegin; 5015 ierr = PetscFree(*(char**)filenametemplate);CHKERRQ(ierr); 5016 PetscFunctionReturn(0); 5017 } 5018 5019 #undef __FUNCT__ 5020 #define __FUNCT__ "TSGetAdapt" 5021 /*@ 5022 TSGetAdapt - Get the adaptive controller context for the current method 5023 5024 Collective on TS if controller has not been created yet 5025 5026 Input Arguments: 5027 . ts - time stepping context 5028 5029 Output Arguments: 5030 . adapt - adaptive controller 5031 5032 Level: intermediate 5033 5034 .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose() 5035 @*/ 5036 PetscErrorCode TSGetAdapt(TS ts,TSAdapt *adapt) 5037 { 5038 PetscErrorCode ierr; 5039 5040 PetscFunctionBegin; 5041 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5042 PetscValidPointer(adapt,2); 5043 if (!ts->adapt) { 5044 ierr = TSAdaptCreate(PetscObjectComm((PetscObject)ts),&ts->adapt);CHKERRQ(ierr); 5045 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->adapt);CHKERRQ(ierr); 5046 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);CHKERRQ(ierr); 5047 } 5048 *adapt = ts->adapt; 5049 PetscFunctionReturn(0); 5050 } 5051 5052 #undef __FUNCT__ 5053 #define __FUNCT__ "TSSetTolerances" 5054 /*@ 5055 TSSetTolerances - Set tolerances for local truncation error when using adaptive controller 5056 5057 Logically Collective 5058 5059 Input Arguments: 5060 + ts - time integration context 5061 . atol - scalar absolute tolerances, PETSC_DECIDE to leave current value 5062 . vatol - vector of absolute tolerances or NULL, used in preference to atol if present 5063 . rtol - scalar relative tolerances, PETSC_DECIDE to leave current value 5064 - vrtol - vector of relative tolerances or NULL, used in preference to atol if present 5065 5066 Options Database keys: 5067 + -ts_rtol <rtol> - relative tolerance for local truncation error 5068 - -ts_atol <atol> Absolute tolerance for local truncation error 5069 5070 Notes: 5071 With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error 5072 (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be 5073 computed only for the differential or the algebraic part then this can be done using the vector of 5074 tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the 5075 differential part and infinity for the algebraic part, the LTE calculation will include only the 5076 differential variables. 5077 5078 Level: beginner 5079 5080 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSGetTolerances() 5081 @*/ 5082 PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol) 5083 { 5084 PetscErrorCode ierr; 5085 5086 PetscFunctionBegin; 5087 if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol; 5088 if (vatol) { 5089 ierr = PetscObjectReference((PetscObject)vatol);CHKERRQ(ierr); 5090 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 5091 5092 ts->vatol = vatol; 5093 } 5094 if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol; 5095 if (vrtol) { 5096 ierr = PetscObjectReference((PetscObject)vrtol);CHKERRQ(ierr); 5097 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 5098 5099 ts->vrtol = vrtol; 5100 } 5101 PetscFunctionReturn(0); 5102 } 5103 5104 #undef __FUNCT__ 5105 #define __FUNCT__ "TSGetTolerances" 5106 /*@ 5107 TSGetTolerances - Get tolerances for local truncation error when using adaptive controller 5108 5109 Logically Collective 5110 5111 Input Arguments: 5112 . ts - time integration context 5113 5114 Output Arguments: 5115 + atol - scalar absolute tolerances, NULL to ignore 5116 . vatol - vector of absolute tolerances, NULL to ignore 5117 . rtol - scalar relative tolerances, NULL to ignore 5118 - vrtol - vector of relative tolerances, NULL to ignore 5119 5120 Level: beginner 5121 5122 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSSetTolerances() 5123 @*/ 5124 PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol) 5125 { 5126 PetscFunctionBegin; 5127 if (atol) *atol = ts->atol; 5128 if (vatol) *vatol = ts->vatol; 5129 if (rtol) *rtol = ts->rtol; 5130 if (vrtol) *vrtol = ts->vrtol; 5131 PetscFunctionReturn(0); 5132 } 5133 5134 #undef __FUNCT__ 5135 #define __FUNCT__ "TSErrorWeightedNorm2" 5136 /*@ 5137 TSErrorWeightedNorm2 - compute a weighted 2-norm of the difference between two state vectors 5138 5139 Collective on TS 5140 5141 Input Arguments: 5142 + ts - time stepping context 5143 . U - state vector, usually ts->vec_sol 5144 - Y - state vector to be compared to U 5145 5146 Output Arguments: 5147 . norm - weighted norm, a value of 1.0 is considered small 5148 5149 Level: developer 5150 5151 .seealso: TSErrorWeightedNorm(), TSErrorWeightedNormInfinity() 5152 @*/ 5153 PetscErrorCode TSErrorWeightedNorm2(TS ts,Vec U,Vec Y,PetscReal *norm) 5154 { 5155 PetscErrorCode ierr; 5156 PetscInt i,n,N,rstart; 5157 const PetscScalar *u,*y; 5158 PetscReal sum,gsum; 5159 PetscReal tol; 5160 5161 PetscFunctionBegin; 5162 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5163 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 5164 PetscValidHeaderSpecific(Y,VEC_CLASSID,3); 5165 PetscValidType(U,2); 5166 PetscValidType(Y,3); 5167 PetscCheckSameComm(U,2,Y,3); 5168 PetscValidPointer(norm,4); 5169 if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"U and Y cannot be the same vector"); 5170 5171 ierr = VecGetSize(U,&N);CHKERRQ(ierr); 5172 ierr = VecGetLocalSize(U,&n);CHKERRQ(ierr); 5173 ierr = VecGetOwnershipRange(U,&rstart,NULL);CHKERRQ(ierr); 5174 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 5175 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 5176 sum = 0.; 5177 if (ts->vatol && ts->vrtol) { 5178 const PetscScalar *atol,*rtol; 5179 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 5180 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 5181 for (i=0; i<n; i++) { 5182 tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 5183 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 5184 } 5185 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 5186 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 5187 } else if (ts->vatol) { /* vector atol, scalar rtol */ 5188 const PetscScalar *atol; 5189 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 5190 for (i=0; i<n; i++) { 5191 tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 5192 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 5193 } 5194 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 5195 } else if (ts->vrtol) { /* scalar atol, vector rtol */ 5196 const PetscScalar *rtol; 5197 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 5198 for (i=0; i<n; i++) { 5199 tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 5200 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 5201 } 5202 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 5203 } else { /* scalar atol, scalar rtol */ 5204 for (i=0; i<n; i++) { 5205 tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 5206 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 5207 } 5208 } 5209 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 5210 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 5211 5212 ierr = MPI_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 5213 *norm = PetscSqrtReal(gsum / N); 5214 5215 if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 5216 PetscFunctionReturn(0); 5217 } 5218 5219 #undef __FUNCT__ 5220 #define __FUNCT__ "TSErrorWeightedNormInfinity" 5221 /*@ 5222 TSErrorWeightedNormInfinity - compute a weighted infinity-norm of the difference between two state vectors 5223 5224 Collective on TS 5225 5226 Input Arguments: 5227 + ts - time stepping context 5228 . U - state vector, usually ts->vec_sol 5229 - Y - state vector to be compared to U 5230 5231 Output Arguments: 5232 . norm - weighted norm, a value of 1.0 is considered small 5233 5234 Level: developer 5235 5236 .seealso: TSErrorWeightedNorm(), TSErrorWeightedNorm2() 5237 @*/ 5238 PetscErrorCode TSErrorWeightedNormInfinity(TS ts,Vec U,Vec Y,PetscReal *norm) 5239 { 5240 PetscErrorCode ierr; 5241 PetscInt i,n,N,rstart,k; 5242 const PetscScalar *u,*y; 5243 PetscReal max,gmax; 5244 PetscReal tol; 5245 5246 PetscFunctionBegin; 5247 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5248 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 5249 PetscValidHeaderSpecific(Y,VEC_CLASSID,3); 5250 PetscValidType(U,2); 5251 PetscValidType(Y,3); 5252 PetscCheckSameComm(U,2,Y,3); 5253 PetscValidPointer(norm,4); 5254 if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"U and Y cannot be the same vector"); 5255 5256 ierr = VecGetSize(U,&N);CHKERRQ(ierr); 5257 ierr = VecGetLocalSize(U,&n);CHKERRQ(ierr); 5258 ierr = VecGetOwnershipRange(U,&rstart,NULL);CHKERRQ(ierr); 5259 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 5260 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 5261 if (ts->vatol && ts->vrtol) { 5262 const PetscScalar *atol,*rtol; 5263 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 5264 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 5265 k = 0; 5266 tol = PetscRealPart(atol[k]) + PetscRealPart(rtol[k]) * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5267 max = PetscAbsScalar(y[k] - u[k]) / tol; 5268 for (i=1; i<n; i++) { 5269 tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 5270 max = PetscMax(max,PetscAbsScalar(y[i] - u[i]) / tol); 5271 } 5272 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 5273 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 5274 } else if (ts->vatol) { /* vector atol, scalar rtol */ 5275 const PetscScalar *atol; 5276 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 5277 k = 0; 5278 tol = PetscRealPart(atol[k]) + ts->rtol * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5279 max = PetscAbsScalar(y[k] - u[k]) / tol; 5280 for (i=1; i<n; i++) { 5281 tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 5282 max = PetscMax(max,PetscAbsScalar(y[i] - u[i]) / tol); 5283 } 5284 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 5285 } else if (ts->vrtol) { /* scalar atol, vector rtol */ 5286 const PetscScalar *rtol; 5287 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 5288 k = 0; 5289 tol = ts->atol + PetscRealPart(rtol[k]) * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5290 max = PetscAbsScalar(y[k] - u[k]) / tol; 5291 for (i=1; i<n; i++) { 5292 tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 5293 max = PetscMax(max,PetscAbsScalar(y[i] - u[i]) / tol); 5294 } 5295 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 5296 } else { /* scalar atol, scalar rtol */ 5297 k = 0; 5298 tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k])); 5299 max = PetscAbsScalar(y[k] - u[k]) / tol; 5300 for (i=1; i<n; i++) { 5301 tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 5302 max = PetscMax(max,PetscAbsScalar(y[i] - u[i]) / tol); 5303 } 5304 } 5305 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 5306 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 5307 5308 ierr = MPI_Allreduce(&max,&gmax,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 5309 *norm = gmax; 5310 5311 if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 5312 PetscFunctionReturn(0); 5313 } 5314 5315 #undef __FUNCT__ 5316 #define __FUNCT__ "TSErrorWeightedNorm" 5317 /*@ 5318 TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors 5319 5320 Collective on TS 5321 5322 Input Arguments: 5323 + ts - time stepping context 5324 . U - state vector, usually ts->vec_sol 5325 . Y - state vector to be compared to U 5326 - wnormtype - norm type, either NORM_2 or NORM_INFINITY 5327 5328 Output Arguments: 5329 . norm - weighted norm, a value of 1.0 is considered small 5330 5331 5332 Options Database Keys: 5333 . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY 5334 5335 Level: developer 5336 5337 .seealso: TSErrorWeightedNormInfinity(), TSErrorWeightedNorm2() 5338 @*/ 5339 PetscErrorCode TSErrorWeightedNorm(TS ts,Vec U,Vec Y,NormType wnormtype,PetscReal *norm) 5340 { 5341 PetscErrorCode ierr; 5342 5343 PetscFunctionBegin; 5344 if (wnormtype == NORM_2) { 5345 ierr = TSErrorWeightedNorm2(ts,U,Y,norm);CHKERRQ(ierr); 5346 } else if(wnormtype == NORM_INFINITY) { 5347 ierr = TSErrorWeightedNormInfinity(ts,U,Y,norm);CHKERRQ(ierr); 5348 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for norm type %s",NormTypes[wnormtype]); 5349 PetscFunctionReturn(0); 5350 } 5351 5352 #undef __FUNCT__ 5353 #define __FUNCT__ "TSSetCFLTimeLocal" 5354 /*@ 5355 TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler 5356 5357 Logically Collective on TS 5358 5359 Input Arguments: 5360 + ts - time stepping context 5361 - cfltime - maximum stable time step if using forward Euler (value can be different on each process) 5362 5363 Note: 5364 After calling this function, the global CFL time can be obtained by calling TSGetCFLTime() 5365 5366 Level: intermediate 5367 5368 .seealso: TSGetCFLTime(), TSADAPTCFL 5369 @*/ 5370 PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime) 5371 { 5372 PetscFunctionBegin; 5373 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5374 ts->cfltime_local = cfltime; 5375 ts->cfltime = -1.; 5376 PetscFunctionReturn(0); 5377 } 5378 5379 #undef __FUNCT__ 5380 #define __FUNCT__ "TSGetCFLTime" 5381 /*@ 5382 TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler 5383 5384 Collective on TS 5385 5386 Input Arguments: 5387 . ts - time stepping context 5388 5389 Output Arguments: 5390 . cfltime - maximum stable time step for forward Euler 5391 5392 Level: advanced 5393 5394 .seealso: TSSetCFLTimeLocal() 5395 @*/ 5396 PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime) 5397 { 5398 PetscErrorCode ierr; 5399 5400 PetscFunctionBegin; 5401 if (ts->cfltime < 0) { 5402 ierr = MPI_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 5403 } 5404 *cfltime = ts->cfltime; 5405 PetscFunctionReturn(0); 5406 } 5407 5408 #undef __FUNCT__ 5409 #define __FUNCT__ "TSVISetVariableBounds" 5410 /*@ 5411 TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu 5412 5413 Input Parameters: 5414 . ts - the TS context. 5415 . xl - lower bound. 5416 . xu - upper bound. 5417 5418 Notes: 5419 If this routine is not called then the lower and upper bounds are set to 5420 PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp(). 5421 5422 Level: advanced 5423 5424 @*/ 5425 PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu) 5426 { 5427 PetscErrorCode ierr; 5428 SNES snes; 5429 5430 PetscFunctionBegin; 5431 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 5432 ierr = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(ierr); 5433 PetscFunctionReturn(0); 5434 } 5435 5436 #if defined(PETSC_HAVE_MATLAB_ENGINE) 5437 #include <mex.h> 5438 5439 typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext; 5440 5441 #undef __FUNCT__ 5442 #define __FUNCT__ "TSComputeFunction_Matlab" 5443 /* 5444 TSComputeFunction_Matlab - Calls the function that has been set with 5445 TSSetFunctionMatlab(). 5446 5447 Collective on TS 5448 5449 Input Parameters: 5450 + snes - the TS context 5451 - u - input vector 5452 5453 Output Parameter: 5454 . y - function vector, as set by TSSetFunction() 5455 5456 Notes: 5457 TSComputeFunction() is typically used within nonlinear solvers 5458 implementations, so most users would not generally call this routine 5459 themselves. 5460 5461 Level: developer 5462 5463 .keywords: TS, nonlinear, compute, function 5464 5465 .seealso: TSSetFunction(), TSGetFunction() 5466 */ 5467 PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec u,Vec udot,Vec y, void *ctx) 5468 { 5469 PetscErrorCode ierr; 5470 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 5471 int nlhs = 1,nrhs = 7; 5472 mxArray *plhs[1],*prhs[7]; 5473 long long int lx = 0,lxdot = 0,ly = 0,ls = 0; 5474 5475 PetscFunctionBegin; 5476 PetscValidHeaderSpecific(snes,TS_CLASSID,1); 5477 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 5478 PetscValidHeaderSpecific(udot,VEC_CLASSID,4); 5479 PetscValidHeaderSpecific(y,VEC_CLASSID,5); 5480 PetscCheckSameComm(snes,1,u,3); 5481 PetscCheckSameComm(snes,1,y,5); 5482 5483 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5484 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 5485 ierr = PetscMemcpy(&lxdot,&udot,sizeof(udot));CHKERRQ(ierr); 5486 ierr = PetscMemcpy(&ly,&y,sizeof(u));CHKERRQ(ierr); 5487 5488 prhs[0] = mxCreateDoubleScalar((double)ls); 5489 prhs[1] = mxCreateDoubleScalar(time); 5490 prhs[2] = mxCreateDoubleScalar((double)lx); 5491 prhs[3] = mxCreateDoubleScalar((double)lxdot); 5492 prhs[4] = mxCreateDoubleScalar((double)ly); 5493 prhs[5] = mxCreateString(sctx->funcname); 5494 prhs[6] = sctx->ctx; 5495 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");CHKERRQ(ierr); 5496 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5497 mxDestroyArray(prhs[0]); 5498 mxDestroyArray(prhs[1]); 5499 mxDestroyArray(prhs[2]); 5500 mxDestroyArray(prhs[3]); 5501 mxDestroyArray(prhs[4]); 5502 mxDestroyArray(prhs[5]); 5503 mxDestroyArray(plhs[0]); 5504 PetscFunctionReturn(0); 5505 } 5506 5507 5508 #undef __FUNCT__ 5509 #define __FUNCT__ "TSSetFunctionMatlab" 5510 /* 5511 TSSetFunctionMatlab - Sets the function evaluation routine and function 5512 vector for use by the TS routines in solving ODEs 5513 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 5514 5515 Logically Collective on TS 5516 5517 Input Parameters: 5518 + ts - the TS context 5519 - func - function evaluation routine 5520 5521 Calling sequence of func: 5522 $ func (TS ts,PetscReal time,Vec u,Vec udot,Vec f,void *ctx); 5523 5524 Level: beginner 5525 5526 .keywords: TS, nonlinear, set, function 5527 5528 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 5529 */ 5530 PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx) 5531 { 5532 PetscErrorCode ierr; 5533 TSMatlabContext *sctx; 5534 5535 PetscFunctionBegin; 5536 /* currently sctx is memory bleed */ 5537 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 5538 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5539 /* 5540 This should work, but it doesn't 5541 sctx->ctx = ctx; 5542 mexMakeArrayPersistent(sctx->ctx); 5543 */ 5544 sctx->ctx = mxDuplicateArray(ctx); 5545 5546 ierr = TSSetIFunction(ts,NULL,TSComputeFunction_Matlab,sctx);CHKERRQ(ierr); 5547 PetscFunctionReturn(0); 5548 } 5549 5550 #undef __FUNCT__ 5551 #define __FUNCT__ "TSComputeJacobian_Matlab" 5552 /* 5553 TSComputeJacobian_Matlab - Calls the function that has been set with 5554 TSSetJacobianMatlab(). 5555 5556 Collective on TS 5557 5558 Input Parameters: 5559 + ts - the TS context 5560 . u - input vector 5561 . A, B - the matrices 5562 - ctx - user context 5563 5564 Level: developer 5565 5566 .keywords: TS, nonlinear, compute, function 5567 5568 .seealso: TSSetFunction(), TSGetFunction() 5569 @*/ 5570 PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec u,Vec udot,PetscReal shift,Mat A,Mat B,void *ctx) 5571 { 5572 PetscErrorCode ierr; 5573 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 5574 int nlhs = 2,nrhs = 9; 5575 mxArray *plhs[2],*prhs[9]; 5576 long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0; 5577 5578 PetscFunctionBegin; 5579 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5580 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 5581 5582 /* call Matlab function in ctx with arguments u and y */ 5583 5584 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 5585 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 5586 ierr = PetscMemcpy(&lxdot,&udot,sizeof(u));CHKERRQ(ierr); 5587 ierr = PetscMemcpy(&lA,A,sizeof(u));CHKERRQ(ierr); 5588 ierr = PetscMemcpy(&lB,B,sizeof(u));CHKERRQ(ierr); 5589 5590 prhs[0] = mxCreateDoubleScalar((double)ls); 5591 prhs[1] = mxCreateDoubleScalar((double)time); 5592 prhs[2] = mxCreateDoubleScalar((double)lx); 5593 prhs[3] = mxCreateDoubleScalar((double)lxdot); 5594 prhs[4] = mxCreateDoubleScalar((double)shift); 5595 prhs[5] = mxCreateDoubleScalar((double)lA); 5596 prhs[6] = mxCreateDoubleScalar((double)lB); 5597 prhs[7] = mxCreateString(sctx->funcname); 5598 prhs[8] = sctx->ctx; 5599 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");CHKERRQ(ierr); 5600 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5601 mxDestroyArray(prhs[0]); 5602 mxDestroyArray(prhs[1]); 5603 mxDestroyArray(prhs[2]); 5604 mxDestroyArray(prhs[3]); 5605 mxDestroyArray(prhs[4]); 5606 mxDestroyArray(prhs[5]); 5607 mxDestroyArray(prhs[6]); 5608 mxDestroyArray(prhs[7]); 5609 mxDestroyArray(plhs[0]); 5610 mxDestroyArray(plhs[1]); 5611 PetscFunctionReturn(0); 5612 } 5613 5614 5615 #undef __FUNCT__ 5616 #define __FUNCT__ "TSSetJacobianMatlab" 5617 /* 5618 TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 5619 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 5620 5621 Logically Collective on TS 5622 5623 Input Parameters: 5624 + ts - the TS context 5625 . A,B - Jacobian matrices 5626 . func - function evaluation routine 5627 - ctx - user context 5628 5629 Calling sequence of func: 5630 $ flag = func (TS ts,PetscReal time,Vec u,Vec udot,Mat A,Mat B,void *ctx); 5631 5632 5633 Level: developer 5634 5635 .keywords: TS, nonlinear, set, function 5636 5637 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 5638 */ 5639 PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx) 5640 { 5641 PetscErrorCode ierr; 5642 TSMatlabContext *sctx; 5643 5644 PetscFunctionBegin; 5645 /* currently sctx is memory bleed */ 5646 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 5647 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5648 /* 5649 This should work, but it doesn't 5650 sctx->ctx = ctx; 5651 mexMakeArrayPersistent(sctx->ctx); 5652 */ 5653 sctx->ctx = mxDuplicateArray(ctx); 5654 5655 ierr = TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 5656 PetscFunctionReturn(0); 5657 } 5658 5659 #undef __FUNCT__ 5660 #define __FUNCT__ "TSMonitor_Matlab" 5661 /* 5662 TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab(). 5663 5664 Collective on TS 5665 5666 .seealso: TSSetFunction(), TSGetFunction() 5667 @*/ 5668 PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec u, void *ctx) 5669 { 5670 PetscErrorCode ierr; 5671 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 5672 int nlhs = 1,nrhs = 6; 5673 mxArray *plhs[1],*prhs[6]; 5674 long long int lx = 0,ls = 0; 5675 5676 PetscFunctionBegin; 5677 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5678 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 5679 5680 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 5681 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 5682 5683 prhs[0] = mxCreateDoubleScalar((double)ls); 5684 prhs[1] = mxCreateDoubleScalar((double)it); 5685 prhs[2] = mxCreateDoubleScalar((double)time); 5686 prhs[3] = mxCreateDoubleScalar((double)lx); 5687 prhs[4] = mxCreateString(sctx->funcname); 5688 prhs[5] = sctx->ctx; 5689 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");CHKERRQ(ierr); 5690 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5691 mxDestroyArray(prhs[0]); 5692 mxDestroyArray(prhs[1]); 5693 mxDestroyArray(prhs[2]); 5694 mxDestroyArray(prhs[3]); 5695 mxDestroyArray(prhs[4]); 5696 mxDestroyArray(plhs[0]); 5697 PetscFunctionReturn(0); 5698 } 5699 5700 5701 #undef __FUNCT__ 5702 #define __FUNCT__ "TSMonitorSetMatlab" 5703 /* 5704 TSMonitorSetMatlab - Sets the monitor function from Matlab 5705 5706 Level: developer 5707 5708 .keywords: TS, nonlinear, set, function 5709 5710 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 5711 */ 5712 PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx) 5713 { 5714 PetscErrorCode ierr; 5715 TSMatlabContext *sctx; 5716 5717 PetscFunctionBegin; 5718 /* currently sctx is memory bleed */ 5719 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 5720 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5721 /* 5722 This should work, but it doesn't 5723 sctx->ctx = ctx; 5724 mexMakeArrayPersistent(sctx->ctx); 5725 */ 5726 sctx->ctx = mxDuplicateArray(ctx); 5727 5728 ierr = TSMonitorSet(ts,TSMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 5729 PetscFunctionReturn(0); 5730 } 5731 #endif 5732 5733 #undef __FUNCT__ 5734 #define __FUNCT__ "TSMonitorLGSolution" 5735 /*@C 5736 TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector 5737 in a time based line graph 5738 5739 Collective on TS 5740 5741 Input Parameters: 5742 + ts - the TS context 5743 . step - current time-step 5744 . ptime - current time 5745 - lg - a line graph object 5746 5747 Options Database: 5748 . -ts_monitor_lg_solution_variables 5749 5750 Level: intermediate 5751 5752 Notes: each process in a parallel run displays its component solutions in a separate window 5753 5754 .keywords: TS, vector, monitor, view 5755 5756 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 5757 @*/ 5758 PetscErrorCode TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 5759 { 5760 PetscErrorCode ierr; 5761 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 5762 const PetscScalar *yy; 5763 PetscInt dim; 5764 Vec v; 5765 5766 PetscFunctionBegin; 5767 if (!step) { 5768 PetscDrawAxis axis; 5769 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 5770 ierr = PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");CHKERRQ(ierr); 5771 if (ctx->names && !ctx->displaynames) { 5772 char **displaynames; 5773 PetscBool flg; 5774 5775 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 5776 ierr = PetscMalloc((dim+1)*sizeof(char*),&displaynames);CHKERRQ(ierr); 5777 ierr = PetscMemzero(displaynames,(dim+1)*sizeof(char*));CHKERRQ(ierr); 5778 ierr = PetscOptionsGetStringArray(((PetscObject)ts)->prefix,"-ts_monitor_lg_solution_variables",displaynames,&dim,&flg);CHKERRQ(ierr); 5779 if (flg) { 5780 ierr = TSMonitorLGCtxSetDisplayVariables(ctx,(const char *const *)displaynames);CHKERRQ(ierr); 5781 } 5782 ierr = PetscStrArrayDestroy(&displaynames);CHKERRQ(ierr); 5783 } 5784 if (ctx->displaynames) { 5785 ierr = PetscDrawLGSetDimension(ctx->lg,ctx->ndisplayvariables);CHKERRQ(ierr); 5786 ierr = PetscDrawLGSetLegend(ctx->lg,(const char *const *)ctx->displaynames);CHKERRQ(ierr); 5787 } else if (ctx->names) { 5788 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 5789 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 5790 ierr = PetscDrawLGSetLegend(ctx->lg,(const char *const *)ctx->names);CHKERRQ(ierr); 5791 } 5792 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 5793 } 5794 if (ctx->transform) { 5795 ierr = (*ctx->transform)(ctx->transformctx,u,&v);CHKERRQ(ierr); 5796 } else { 5797 v = u; 5798 } 5799 ierr = VecGetArrayRead(v,&yy);CHKERRQ(ierr); 5800 #if defined(PETSC_USE_COMPLEX) 5801 { 5802 PetscReal *yreal; 5803 PetscInt i,n; 5804 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 5805 ierr = PetscMalloc1(n,&yreal);CHKERRQ(ierr); 5806 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 5807 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 5808 ierr = PetscFree(yreal);CHKERRQ(ierr); 5809 } 5810 #else 5811 if (ctx->displaynames) { 5812 PetscInt i; 5813 for (i=0; i<ctx->ndisplayvariables; i++) { 5814 ctx->displayvalues[i] = yy[ctx->displayvariables[i]]; 5815 } 5816 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,ctx->displayvalues);CHKERRQ(ierr); 5817 } else { 5818 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 5819 } 5820 #endif 5821 ierr = VecRestoreArrayRead(v,&yy);CHKERRQ(ierr); 5822 if (ctx->transform) { 5823 ierr = VecDestroy(&v);CHKERRQ(ierr); 5824 } 5825 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 5826 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 5827 } 5828 PetscFunctionReturn(0); 5829 } 5830 5831 5832 #undef __FUNCT__ 5833 #define __FUNCT__ "TSMonitorLGSetVariableNames" 5834 /*@C 5835 TSMonitorLGSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot 5836 5837 Collective on TS 5838 5839 Input Parameters: 5840 + ts - the TS context 5841 - names - the names of the components, final string must be NULL 5842 5843 Level: intermediate 5844 5845 .keywords: TS, vector, monitor, view 5846 5847 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetVariableNames() 5848 @*/ 5849 PetscErrorCode TSMonitorLGSetVariableNames(TS ts,const char * const *names) 5850 { 5851 PetscErrorCode ierr; 5852 PetscInt i; 5853 5854 PetscFunctionBegin; 5855 for (i=0; i<ts->numbermonitors; i++) { 5856 if (ts->monitor[i] == TSMonitorLGSolution) { 5857 ierr = TSMonitorLGCtxSetVariableNames((TSMonitorLGCtx)ts->monitorcontext[i],names);CHKERRQ(ierr); 5858 break; 5859 } 5860 } 5861 PetscFunctionReturn(0); 5862 } 5863 5864 #undef __FUNCT__ 5865 #define __FUNCT__ "TSMonitorLGCtxSetVariableNames" 5866 /*@C 5867 TSMonitorLGCtxSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot 5868 5869 Collective on TS 5870 5871 Input Parameters: 5872 + ts - the TS context 5873 - names - the names of the components, final string must be NULL 5874 5875 Level: intermediate 5876 5877 .keywords: TS, vector, monitor, view 5878 5879 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables(), TSMonitorLGSetVariableNames() 5880 @*/ 5881 PetscErrorCode TSMonitorLGCtxSetVariableNames(TSMonitorLGCtx ctx,const char * const *names) 5882 { 5883 PetscErrorCode ierr; 5884 5885 PetscFunctionBegin; 5886 ierr = PetscStrArrayDestroy(&ctx->names);CHKERRQ(ierr); 5887 ierr = PetscStrArrayallocpy(names,&ctx->names);CHKERRQ(ierr); 5888 PetscFunctionReturn(0); 5889 } 5890 5891 #undef __FUNCT__ 5892 #define __FUNCT__ "TSMonitorLGGetVariableNames" 5893 /*@C 5894 TSMonitorLGGetVariableNames - Gets the name of each component in the solution vector so that it may be displayed in the plot 5895 5896 Collective on TS 5897 5898 Input Parameter: 5899 . ts - the TS context 5900 5901 Output Parameter: 5902 . names - the names of the components, final string must be NULL 5903 5904 Level: intermediate 5905 5906 .keywords: TS, vector, monitor, view 5907 5908 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables() 5909 @*/ 5910 PetscErrorCode TSMonitorLGGetVariableNames(TS ts,const char *const **names) 5911 { 5912 PetscInt i; 5913 5914 PetscFunctionBegin; 5915 *names = NULL; 5916 for (i=0; i<ts->numbermonitors; i++) { 5917 if (ts->monitor[i] == TSMonitorLGSolution) { 5918 TSMonitorLGCtx ctx = (TSMonitorLGCtx) ts->monitorcontext[i]; 5919 *names = (const char *const *)ctx->names; 5920 break; 5921 } 5922 } 5923 PetscFunctionReturn(0); 5924 } 5925 5926 #undef __FUNCT__ 5927 #define __FUNCT__ "TSMonitorLGCtxSetDisplayVariables" 5928 /*@C 5929 TSMonitorLGCtxSetDisplayVariables - Sets the variables that are to be display in the monitor 5930 5931 Collective on TS 5932 5933 Input Parameters: 5934 + ctx - the TSMonitorLG context 5935 . displaynames - the names of the components, final string must be NULL 5936 5937 Level: intermediate 5938 5939 .keywords: TS, vector, monitor, view 5940 5941 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames() 5942 @*/ 5943 PetscErrorCode TSMonitorLGCtxSetDisplayVariables(TSMonitorLGCtx ctx,const char * const *displaynames) 5944 { 5945 PetscInt j = 0,k; 5946 PetscErrorCode ierr; 5947 5948 PetscFunctionBegin; 5949 if (!ctx->names) PetscFunctionReturn(0); 5950 ierr = PetscStrArrayDestroy(&ctx->displaynames);CHKERRQ(ierr); 5951 ierr = PetscStrArrayallocpy(displaynames,&ctx->displaynames);CHKERRQ(ierr); 5952 while (displaynames[j]) j++; 5953 ctx->ndisplayvariables = j; 5954 ierr = PetscMalloc1(ctx->ndisplayvariables,&ctx->displayvariables);CHKERRQ(ierr); 5955 ierr = PetscMalloc1(ctx->ndisplayvariables,&ctx->displayvalues);CHKERRQ(ierr); 5956 j = 0; 5957 while (displaynames[j]) { 5958 k = 0; 5959 while (ctx->names[k]) { 5960 PetscBool flg; 5961 ierr = PetscStrcmp(displaynames[j],ctx->names[k],&flg);CHKERRQ(ierr); 5962 if (flg) { 5963 ctx->displayvariables[j] = k; 5964 break; 5965 } 5966 k++; 5967 } 5968 j++; 5969 } 5970 PetscFunctionReturn(0); 5971 } 5972 5973 5974 #undef __FUNCT__ 5975 #define __FUNCT__ "TSMonitorLGSetDisplayVariables" 5976 /*@C 5977 TSMonitorLGSetDisplayVariables - Sets the variables that are to be display in the monitor 5978 5979 Collective on TS 5980 5981 Input Parameters: 5982 + ts - the TS context 5983 . displaynames - the names of the components, final string must be NULL 5984 5985 Level: intermediate 5986 5987 .keywords: TS, vector, monitor, view 5988 5989 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames() 5990 @*/ 5991 PetscErrorCode TSMonitorLGSetDisplayVariables(TS ts,const char * const *displaynames) 5992 { 5993 PetscInt i; 5994 PetscErrorCode ierr; 5995 5996 PetscFunctionBegin; 5997 for (i=0; i<ts->numbermonitors; i++) { 5998 if (ts->monitor[i] == TSMonitorLGSolution) { 5999 ierr = TSMonitorLGCtxSetDisplayVariables((TSMonitorLGCtx)ts->monitorcontext[i],displaynames);CHKERRQ(ierr); 6000 break; 6001 } 6002 } 6003 PetscFunctionReturn(0); 6004 } 6005 6006 #undef __FUNCT__ 6007 #define __FUNCT__ "TSMonitorLGSetTransform" 6008 /*@C 6009 TSMonitorLGSetTransform - Solution vector will be transformed by provided function before being displayed 6010 6011 Collective on TS 6012 6013 Input Parameters: 6014 + ts - the TS context 6015 . transform - the transform function 6016 . destroy - function to destroy the optional context 6017 - ctx - optional context used by transform function 6018 6019 Level: intermediate 6020 6021 .keywords: TS, vector, monitor, view 6022 6023 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames(), TSMonitorLGCtxSetTransform() 6024 @*/ 6025 PetscErrorCode TSMonitorLGSetTransform(TS ts,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx) 6026 { 6027 PetscInt i; 6028 PetscErrorCode ierr; 6029 6030 PetscFunctionBegin; 6031 for (i=0; i<ts->numbermonitors; i++) { 6032 if (ts->monitor[i] == TSMonitorLGSolution) { 6033 ierr = TSMonitorLGCtxSetTransform((TSMonitorLGCtx)ts->monitorcontext[i],transform,destroy,tctx);CHKERRQ(ierr); 6034 } 6035 } 6036 PetscFunctionReturn(0); 6037 } 6038 6039 #undef __FUNCT__ 6040 #define __FUNCT__ "TSMonitorLGCtxSetTransform" 6041 /*@C 6042 TSMonitorLGCtxSetTransform - Solution vector will be transformed by provided function before being displayed 6043 6044 Collective on TSLGCtx 6045 6046 Input Parameters: 6047 + ts - the TS context 6048 . transform - the transform function 6049 . destroy - function to destroy the optional context 6050 - ctx - optional context used by transform function 6051 6052 Level: intermediate 6053 6054 .keywords: TS, vector, monitor, view 6055 6056 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames(), TSMonitorLGSetTransform() 6057 @*/ 6058 PetscErrorCode TSMonitorLGCtxSetTransform(TSMonitorLGCtx ctx,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx) 6059 { 6060 PetscFunctionBegin; 6061 ctx->transform = transform; 6062 ctx->transformdestroy = destroy; 6063 ctx->transformctx = tctx; 6064 PetscFunctionReturn(0); 6065 } 6066 6067 #undef __FUNCT__ 6068 #define __FUNCT__ "TSMonitorLGError" 6069 /*@C 6070 TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the solution vector 6071 in a time based line graph 6072 6073 Collective on TS 6074 6075 Input Parameters: 6076 + ts - the TS context 6077 . step - current time-step 6078 . ptime - current time 6079 - lg - a line graph object 6080 6081 Level: intermediate 6082 6083 Notes: 6084 Only for sequential solves. 6085 6086 The user must provide the solution using TSSetSolutionFunction() to use this monitor. 6087 6088 Options Database Keys: 6089 . -ts_monitor_lg_error - create a graphical monitor of error history 6090 6091 .keywords: TS, vector, monitor, view 6092 6093 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction() 6094 @*/ 6095 PetscErrorCode TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 6096 { 6097 PetscErrorCode ierr; 6098 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 6099 const PetscScalar *yy; 6100 Vec y; 6101 PetscInt dim; 6102 6103 PetscFunctionBegin; 6104 if (!step) { 6105 PetscDrawAxis axis; 6106 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 6107 ierr = PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Solution");CHKERRQ(ierr); 6108 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 6109 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 6110 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 6111 } 6112 ierr = VecDuplicate(u,&y);CHKERRQ(ierr); 6113 ierr = TSComputeSolutionFunction(ts,ptime,y);CHKERRQ(ierr); 6114 ierr = VecAXPY(y,-1.0,u);CHKERRQ(ierr); 6115 ierr = VecGetArrayRead(y,&yy);CHKERRQ(ierr); 6116 #if defined(PETSC_USE_COMPLEX) 6117 { 6118 PetscReal *yreal; 6119 PetscInt i,n; 6120 ierr = VecGetLocalSize(y,&n);CHKERRQ(ierr); 6121 ierr = PetscMalloc1(n,&yreal);CHKERRQ(ierr); 6122 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 6123 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 6124 ierr = PetscFree(yreal);CHKERRQ(ierr); 6125 } 6126 #else 6127 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 6128 #endif 6129 ierr = VecRestoreArrayRead(y,&yy);CHKERRQ(ierr); 6130 ierr = VecDestroy(&y);CHKERRQ(ierr); 6131 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 6132 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 6133 } 6134 PetscFunctionReturn(0); 6135 } 6136 6137 #undef __FUNCT__ 6138 #define __FUNCT__ "TSMonitorLGSNESIterations" 6139 PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 6140 { 6141 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 6142 PetscReal x = ptime,y; 6143 PetscErrorCode ierr; 6144 PetscInt its; 6145 6146 PetscFunctionBegin; 6147 if (!n) { 6148 PetscDrawAxis axis; 6149 6150 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 6151 ierr = PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");CHKERRQ(ierr); 6152 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 6153 6154 ctx->snes_its = 0; 6155 } 6156 ierr = TSGetSNESIterations(ts,&its);CHKERRQ(ierr); 6157 y = its - ctx->snes_its; 6158 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 6159 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 6160 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 6161 } 6162 ctx->snes_its = its; 6163 PetscFunctionReturn(0); 6164 } 6165 6166 #undef __FUNCT__ 6167 #define __FUNCT__ "TSMonitorLGKSPIterations" 6168 PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 6169 { 6170 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 6171 PetscReal x = ptime,y; 6172 PetscErrorCode ierr; 6173 PetscInt its; 6174 6175 PetscFunctionBegin; 6176 if (!n) { 6177 PetscDrawAxis axis; 6178 6179 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 6180 ierr = PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");CHKERRQ(ierr); 6181 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 6182 6183 ctx->ksp_its = 0; 6184 } 6185 ierr = TSGetKSPIterations(ts,&its);CHKERRQ(ierr); 6186 y = its - ctx->ksp_its; 6187 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 6188 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 6189 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 6190 } 6191 ctx->ksp_its = its; 6192 PetscFunctionReturn(0); 6193 } 6194 6195 #undef __FUNCT__ 6196 #define __FUNCT__ "TSComputeLinearStability" 6197 /*@ 6198 TSComputeLinearStability - computes the linear stability function at a point 6199 6200 Collective on TS and Vec 6201 6202 Input Parameters: 6203 + ts - the TS context 6204 - xr,xi - real and imaginary part of input arguments 6205 6206 Output Parameters: 6207 . yr,yi - real and imaginary part of function value 6208 6209 Level: developer 6210 6211 .keywords: TS, compute 6212 6213 .seealso: TSSetRHSFunction(), TSComputeIFunction() 6214 @*/ 6215 PetscErrorCode TSComputeLinearStability(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi) 6216 { 6217 PetscErrorCode ierr; 6218 6219 PetscFunctionBegin; 6220 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 6221 if (!ts->ops->linearstability) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Linearized stability function not provided for this method"); 6222 ierr = (*ts->ops->linearstability)(ts,xr,xi,yr,yi);CHKERRQ(ierr); 6223 PetscFunctionReturn(0); 6224 } 6225 6226 /* ------------------------------------------------------------------------*/ 6227 #undef __FUNCT__ 6228 #define __FUNCT__ "TSMonitorEnvelopeCtxCreate" 6229 /*@C 6230 TSMonitorEnvelopeCtxCreate - Creates a context for use with TSMonitorEnvelope() 6231 6232 Collective on TS 6233 6234 Input Parameters: 6235 . ts - the ODE solver object 6236 6237 Output Parameter: 6238 . ctx - the context 6239 6240 Level: intermediate 6241 6242 .keywords: TS, monitor, line graph, residual, seealso 6243 6244 .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError() 6245 6246 @*/ 6247 PetscErrorCode TSMonitorEnvelopeCtxCreate(TS ts,TSMonitorEnvelopeCtx *ctx) 6248 { 6249 PetscErrorCode ierr; 6250 6251 PetscFunctionBegin; 6252 ierr = PetscNew(ctx);CHKERRQ(ierr); 6253 PetscFunctionReturn(0); 6254 } 6255 6256 #undef __FUNCT__ 6257 #define __FUNCT__ "TSMonitorEnvelope" 6258 /*@C 6259 TSMonitorEnvelope - Monitors the maximum and minimum value of each component of the solution 6260 6261 Collective on TS 6262 6263 Input Parameters: 6264 + ts - the TS context 6265 . step - current time-step 6266 . ptime - current time 6267 - ctx - the envelope context 6268 6269 Options Database: 6270 . -ts_monitor_envelope 6271 6272 Level: intermediate 6273 6274 Notes: after a solve you can use TSMonitorEnvelopeGetBounds() to access the envelope 6275 6276 .keywords: TS, vector, monitor, view 6277 6278 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorEnvelopeGetBounds() 6279 @*/ 6280 PetscErrorCode TSMonitorEnvelope(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 6281 { 6282 PetscErrorCode ierr; 6283 TSMonitorEnvelopeCtx ctx = (TSMonitorEnvelopeCtx)dummy; 6284 6285 PetscFunctionBegin; 6286 if (!ctx->max) { 6287 ierr = VecDuplicate(u,&ctx->max);CHKERRQ(ierr); 6288 ierr = VecDuplicate(u,&ctx->min);CHKERRQ(ierr); 6289 ierr = VecCopy(u,ctx->max);CHKERRQ(ierr); 6290 ierr = VecCopy(u,ctx->min);CHKERRQ(ierr); 6291 } else { 6292 ierr = VecPointwiseMax(ctx->max,u,ctx->max);CHKERRQ(ierr); 6293 ierr = VecPointwiseMin(ctx->min,u,ctx->min);CHKERRQ(ierr); 6294 } 6295 PetscFunctionReturn(0); 6296 } 6297 6298 6299 #undef __FUNCT__ 6300 #define __FUNCT__ "TSMonitorEnvelopeGetBounds" 6301 /*@C 6302 TSMonitorEnvelopeGetBounds - Gets the bounds for the components of the solution 6303 6304 Collective on TS 6305 6306 Input Parameter: 6307 . ts - the TS context 6308 6309 Output Parameter: 6310 + max - the maximum values 6311 - min - the minimum values 6312 6313 Level: intermediate 6314 6315 .keywords: TS, vector, monitor, view 6316 6317 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables() 6318 @*/ 6319 PetscErrorCode TSMonitorEnvelopeGetBounds(TS ts,Vec *max,Vec *min) 6320 { 6321 PetscInt i; 6322 6323 PetscFunctionBegin; 6324 if (max) *max = NULL; 6325 if (min) *min = NULL; 6326 for (i=0; i<ts->numbermonitors; i++) { 6327 if (ts->monitor[i] == TSMonitorEnvelope) { 6328 TSMonitorEnvelopeCtx ctx = (TSMonitorEnvelopeCtx) ts->monitorcontext[i]; 6329 if (max) *max = ctx->max; 6330 if (min) *min = ctx->min; 6331 break; 6332 } 6333 } 6334 PetscFunctionReturn(0); 6335 } 6336 6337 #undef __FUNCT__ 6338 #define __FUNCT__ "TSMonitorEnvelopeCtxDestroy" 6339 /*@C 6340 TSMonitorEnvelopeCtxDestroy - Destroys a context that was created with TSMonitorEnvelopeCtxCreate(). 6341 6342 Collective on TSMonitorEnvelopeCtx 6343 6344 Input Parameter: 6345 . ctx - the monitor context 6346 6347 Level: intermediate 6348 6349 .keywords: TS, monitor, line graph, destroy 6350 6351 .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep(); 6352 @*/ 6353 PetscErrorCode TSMonitorEnvelopeCtxDestroy(TSMonitorEnvelopeCtx *ctx) 6354 { 6355 PetscErrorCode ierr; 6356 6357 PetscFunctionBegin; 6358 ierr = VecDestroy(&(*ctx)->min);CHKERRQ(ierr); 6359 ierr = VecDestroy(&(*ctx)->max);CHKERRQ(ierr); 6360 ierr = PetscFree(*ctx);CHKERRQ(ierr); 6361 PetscFunctionReturn(0); 6362 } 6363 6364 #undef __FUNCT__ 6365 #define __FUNCT__ "TSRollBack" 6366 /*@ 6367 TSRollBack - Rolls back one time step 6368 6369 Collective on TS 6370 6371 Input Parameter: 6372 . ts - the TS context obtained from TSCreate() 6373 6374 Level: advanced 6375 6376 .keywords: TS, timestep, rollback 6377 6378 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSInterpolate() 6379 @*/ 6380 PetscErrorCode TSRollBack(TS ts) 6381 { 6382 PetscErrorCode ierr; 6383 6384 PetscFunctionBegin; 6385 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 6386 6387 if (!ts->ops->rollback) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRollBack not implemented for type '%s'",((PetscObject)ts)->type_name); 6388 ierr = (*ts->ops->rollback)(ts);CHKERRQ(ierr); 6389 ts->time_step = ts->ptime - ts->ptime_prev; 6390 ts->ptime = ts->ptime_prev; 6391 ts->steprollback = PETSC_TRUE; /* Flag to indicate that the step is rollbacked */ 6392 PetscFunctionReturn(0); 6393 } 6394 6395 #undef __FUNCT__ 6396 #define __FUNCT__ "TSGetStages" 6397 /*@ 6398 TSGetStages - Get the number of stages and stage values 6399 6400 Input Parameter: 6401 . ts - the TS context obtained from TSCreate() 6402 6403 Level: advanced 6404 6405 .keywords: TS, getstages 6406 6407 .seealso: TSCreate() 6408 @*/ 6409 PetscErrorCode TSGetStages(TS ts,PetscInt *ns, Vec **Y) 6410 { 6411 PetscErrorCode ierr; 6412 6413 PetscFunctionBegin; 6414 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 6415 PetscValidPointer(ns,2); 6416 6417 if (!ts->ops->getstages) *ns=0; 6418 else { 6419 ierr = (*ts->ops->getstages)(ts,ns,Y);CHKERRQ(ierr); 6420 } 6421 PetscFunctionReturn(0); 6422 } 6423 6424 #undef __FUNCT__ 6425 #define __FUNCT__ "TSComputeIJacobianDefaultColor" 6426 /*@C 6427 TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity. 6428 6429 Collective on SNES 6430 6431 Input Parameters: 6432 + ts - the TS context 6433 . t - current timestep 6434 . U - state vector 6435 . Udot - time derivative of state vector 6436 . shift - shift to apply, see note below 6437 - ctx - an optional user context 6438 6439 Output Parameters: 6440 + J - Jacobian matrix (not altered in this routine) 6441 - B - newly computed Jacobian matrix to use with preconditioner (generally the same as J) 6442 6443 Level: intermediate 6444 6445 Notes: 6446 If F(t,U,Udot)=0 is the DAE, the required Jacobian is 6447 6448 dF/dU + shift*dF/dUdot 6449 6450 Most users should not need to explicitly call this routine, as it 6451 is used internally within the nonlinear solvers. 6452 6453 This will first try to get the coloring from the DM. If the DM type has no coloring 6454 routine, then it will try to get the coloring from the matrix. This requires that the 6455 matrix have nonzero entries precomputed. 6456 6457 .keywords: TS, finite differences, Jacobian, coloring, sparse 6458 .seealso: TSSetIJacobian(), MatFDColoringCreate(), MatFDColoringSetFunction() 6459 @*/ 6460 PetscErrorCode TSComputeIJacobianDefaultColor(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat J,Mat B,void *ctx) 6461 { 6462 SNES snes; 6463 MatFDColoring color; 6464 PetscBool hascolor, matcolor = PETSC_FALSE; 6465 PetscErrorCode ierr; 6466 6467 PetscFunctionBegin; 6468 ierr = PetscOptionsGetBool(((PetscObject) ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL);CHKERRQ(ierr); 6469 ierr = PetscObjectQuery((PetscObject) B, "TSMatFDColoring", (PetscObject *) &color);CHKERRQ(ierr); 6470 if (!color) { 6471 DM dm; 6472 ISColoring iscoloring; 6473 6474 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 6475 ierr = DMHasColoring(dm, &hascolor);CHKERRQ(ierr); 6476 if (hascolor && !matcolor) { 6477 ierr = DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring);CHKERRQ(ierr); 6478 ierr = MatFDColoringCreate(B, iscoloring, &color);CHKERRQ(ierr); 6479 ierr = MatFDColoringSetFunction(color, (PetscErrorCode (*)(void)) SNESTSFormFunction, (void *) ts);CHKERRQ(ierr); 6480 ierr = MatFDColoringSetFromOptions(color);CHKERRQ(ierr); 6481 ierr = MatFDColoringSetUp(B, iscoloring, color);CHKERRQ(ierr); 6482 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 6483 } else { 6484 MatColoring mc; 6485 6486 ierr = MatColoringCreate(B, &mc);CHKERRQ(ierr); 6487 ierr = MatColoringSetDistance(mc, 2);CHKERRQ(ierr); 6488 ierr = MatColoringSetType(mc, MATCOLORINGSL);CHKERRQ(ierr); 6489 ierr = MatColoringSetFromOptions(mc);CHKERRQ(ierr); 6490 ierr = MatColoringApply(mc, &iscoloring);CHKERRQ(ierr); 6491 ierr = MatColoringDestroy(&mc);CHKERRQ(ierr); 6492 ierr = MatFDColoringCreate(B, iscoloring, &color);CHKERRQ(ierr); 6493 ierr = MatFDColoringSetFunction(color, (PetscErrorCode (*)(void)) SNESTSFormFunction, (void *) ts);CHKERRQ(ierr); 6494 ierr = MatFDColoringSetFromOptions(color);CHKERRQ(ierr); 6495 ierr = MatFDColoringSetUp(B, iscoloring, color);CHKERRQ(ierr); 6496 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 6497 } 6498 ierr = PetscObjectCompose((PetscObject) B, "TSMatFDColoring", (PetscObject) color);CHKERRQ(ierr); 6499 ierr = PetscObjectDereference((PetscObject) color);CHKERRQ(ierr); 6500 } 6501 ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr); 6502 ierr = MatFDColoringApply(B, color, U, snes);CHKERRQ(ierr); 6503 if (J != B) { 6504 ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 6505 ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 6506 } 6507 PetscFunctionReturn(0); 6508 } 6509 6510 #undef __FUNCT__ 6511 #define __FUNCT__ "TSClone" 6512 /*@C 6513 TSClone - This function clones a time step object. 6514 6515 Collective on MPI_Comm 6516 6517 Input Parameter: 6518 . tsin - The input TS 6519 6520 Output Parameter: 6521 . tsout - The output TS (cloned) 6522 6523 Notes: 6524 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. 6525 6526 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); 6527 6528 Level: developer 6529 6530 .keywords: TS, clone 6531 .seealso: TSCreate(), TSSetType(), TSSetUp(), TSDestroy(), TSSetProblemType() 6532 @*/ 6533 PetscErrorCode TSClone(TS tsin, TS *tsout) 6534 { 6535 TS t; 6536 PetscErrorCode ierr; 6537 SNES snes_start; 6538 DM dm; 6539 TSType type; 6540 6541 PetscFunctionBegin; 6542 PetscValidPointer(tsin,1); 6543 *tsout = NULL; 6544 6545 ierr = PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView);CHKERRQ(ierr); 6546 6547 /* General TS description */ 6548 t->numbermonitors = 0; 6549 t->setupcalled = 0; 6550 t->ksp_its = 0; 6551 t->snes_its = 0; 6552 t->nwork = 0; 6553 t->rhsjacobian.time = -1e20; 6554 t->rhsjacobian.scale = 1.; 6555 t->ijacobian.shift = 1.; 6556 6557 ierr = TSGetSNES(tsin,&snes_start); CHKERRQ(ierr); 6558 ierr = TSSetSNES(t,snes_start); CHKERRQ(ierr); 6559 6560 ierr = TSGetDM(tsin,&dm); CHKERRQ(ierr); 6561 ierr = TSSetDM(t,dm); CHKERRQ(ierr); 6562 6563 t->adapt=tsin->adapt; 6564 PetscObjectReference((PetscObject)t->adapt); 6565 6566 t->problem_type = tsin->problem_type; 6567 t->ptime = tsin->ptime; 6568 t->time_step = tsin->time_step; 6569 t->time_step_orig = tsin->time_step_orig; 6570 t->max_time = tsin->max_time; 6571 t->steps = tsin->steps; 6572 t->max_steps = tsin->max_steps; 6573 t->equation_type = tsin->equation_type; 6574 t->atol = tsin->atol; 6575 t->rtol = tsin->rtol; 6576 t->max_snes_failures = tsin->max_snes_failures; 6577 t->max_reject = tsin->max_reject; 6578 t->errorifstepfailed = tsin->errorifstepfailed; 6579 6580 ierr = TSGetType(tsin,&type); CHKERRQ(ierr); 6581 ierr = TSSetType(t,type); CHKERRQ(ierr); 6582 6583 t->vec_sol = NULL; 6584 6585 t->cfltime = tsin->cfltime; 6586 t->cfltime_local = tsin->cfltime_local; 6587 t->exact_final_time = tsin->exact_final_time; 6588 6589 ierr = PetscMemcpy(t->ops,tsin->ops,sizeof(struct _TSOps));CHKERRQ(ierr); 6590 6591 *tsout = t; 6592 PetscFunctionReturn(0); 6593 } 6594