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