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