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