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