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