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