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