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