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