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