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