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 39 . -ts_save_trajectory - checkpoint the solution at each time-step 40 . -ts_max_time <time> - maximum time to compute to 41 . -ts_max_steps <steps> - maximum number of time-steps to take 42 . -ts_init_time <time> - initial time to start computation 43 . -ts_final_time <time> - final time to compute to (deprecated: use -ts_max_time) 44 . -ts_dt <dt> - initial time step 45 . -ts_exact_final_time <stepover,interpolate,matchstep> - whether to stop at the exact given final time and how to compute the solution at that time 46 . -ts_max_snes_failures <maxfailures> - Maximum number of nonlinear solve failures allowed 47 . -ts_max_reject <maxrejects> - Maximum number of step rejections before step fails 48 . -ts_error_if_step_fails <true,false> - Error if no step succeeds 49 . -ts_rtol <rtol> - relative tolerance for local truncation error 50 . -ts_atol <atol> Absolute tolerance for local truncation error 51 . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - test the Jacobian at each iteration against finite difference with RHS function 52 . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - test the Jacobian at each iteration against finite difference with RHS function 53 . -ts_adjoint_solve <yes,no> After solving the ODE/DAE solve the adjoint problem (requires -ts_save_trajectory) 54 . -ts_fd_color - Use finite differences with coloring to compute IJacobian 55 . -ts_monitor - print information at each timestep 56 . -ts_monitor_cancel - Cancel all monitors 57 . -ts_monitor_lg_solution - Monitor solution graphically 58 . -ts_monitor_lg_error - Monitor error graphically 59 . -ts_monitor_error - Monitors norm of error 60 . -ts_monitor_lg_timestep - Monitor timestep size graphically 61 . -ts_monitor_lg_timestep_log - Monitor log timestep size graphically 62 . -ts_monitor_lg_snes_iterations - Monitor number nonlinear iterations for each timestep graphically 63 . -ts_monitor_lg_ksp_iterations - Monitor number nonlinear iterations for each timestep graphically 64 . -ts_monitor_sp_eig - Monitor eigenvalues of linearized operator graphically 65 . -ts_monitor_draw_solution - Monitor solution graphically 66 . -ts_monitor_draw_solution_phase <xleft,yleft,xright,yright> - Monitor solution graphically with phase diagram, requires problem with exactly 2 degrees of freedom 67 . -ts_monitor_draw_error - Monitor error graphically, requires use to have provided TSSetSolutionFunction() 68 . -ts_monitor_solution [ascii binary draw][:filename][:viewerformat] - monitors the solution at each timestep 69 . -ts_monitor_solution_vtk <filename.vts,filename.vtu> - Save each time step to a binary file, use filename-%%03D.vts (filename-%%03D.vtu) 70 - -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time 71 72 Notes: 73 See SNESSetFromOptions() and KSPSetFromOptions() for how to control the nonlinear and linear solves used by the time-stepper. 74 75 Certain SNES options get reset for each new nonlinear solver, for example -snes_lag_jacobian <its> and -snes_lag_preconditioner <its>, in order 76 to retain them over the multiple nonlinear solves that TS uses you mush also provide -snes_lag_jacobian_persists true and 77 -snes_lag_preconditioner_persists true 78 79 Developer Note: 80 We should unify all the -ts_monitor options in the way that -xxx_view has been unified 81 82 Level: beginner 83 84 .seealso: TSGetType() 85 @*/ 86 PetscErrorCode TSSetFromOptions(TS ts) 87 { 88 PetscBool opt,flg,tflg; 89 PetscErrorCode ierr; 90 char monfilename[PETSC_MAX_PATH_LEN]; 91 PetscReal time_step; 92 TSExactFinalTimeOption eftopt; 93 char dir[16]; 94 TSIFunction ifun; 95 const char *defaultType; 96 char typeName[256]; 97 98 PetscFunctionBegin; 99 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 100 101 ierr = TSRegisterAll();CHKERRQ(ierr); 102 ierr = TSGetIFunction(ts,NULL,&ifun,NULL);CHKERRQ(ierr); 103 104 ierr = PetscObjectOptionsBegin((PetscObject)ts);CHKERRQ(ierr); 105 if (((PetscObject)ts)->type_name) defaultType = ((PetscObject)ts)->type_name; 106 else defaultType = ifun ? TSBEULER : TSEULER; 107 ierr = PetscOptionsFList("-ts_type","TS method","TSSetType",TSList,defaultType,typeName,256,&opt);CHKERRQ(ierr); 108 if (opt) { 109 ierr = TSSetType(ts,typeName);CHKERRQ(ierr); 110 } else { 111 ierr = TSSetType(ts,defaultType);CHKERRQ(ierr); 112 } 113 114 /* Handle generic TS options */ 115 ierr = PetscOptionsDeprecated("-ts_final_time","-ts_max_time","3.10",NULL);CHKERRQ(ierr); 116 ierr = PetscOptionsReal("-ts_max_time","Maximum time to run to","TSSetMaxTime",ts->max_time,&ts->max_time,NULL);CHKERRQ(ierr); 117 ierr = PetscOptionsInt("-ts_max_steps","Maximum number of time steps","TSSetMaxSteps",ts->max_steps,&ts->max_steps,NULL);CHKERRQ(ierr); 118 ierr = PetscOptionsReal("-ts_init_time","Initial time","TSSetTime",ts->ptime,&ts->ptime,NULL);CHKERRQ(ierr); 119 ierr = PetscOptionsReal("-ts_dt","Initial time step","TSSetTimeStep",ts->time_step,&time_step,&flg);CHKERRQ(ierr); 120 if (flg) {ierr = TSSetTimeStep(ts,time_step);CHKERRQ(ierr);} 121 ierr = PetscOptionsEnum("-ts_exact_final_time","Option for handling of final time step","TSSetExactFinalTime",TSExactFinalTimeOptions,(PetscEnum)ts->exact_final_time,(PetscEnum*)&eftopt,&flg);CHKERRQ(ierr); 122 if (flg) {ierr = TSSetExactFinalTime(ts,eftopt);CHKERRQ(ierr);} 123 ierr = PetscOptionsInt("-ts_max_snes_failures","Maximum number of nonlinear solve failures","TSSetMaxSNESFailures",ts->max_snes_failures,&ts->max_snes_failures,NULL);CHKERRQ(ierr); 124 ierr = PetscOptionsInt("-ts_max_reject","Maximum number of step rejections before step fails","TSSetMaxStepRejections",ts->max_reject,&ts->max_reject,NULL);CHKERRQ(ierr); 125 ierr = PetscOptionsBool("-ts_error_if_step_fails","Error if no step succeeds","TSSetErrorIfStepFails",ts->errorifstepfailed,&ts->errorifstepfailed,NULL);CHKERRQ(ierr); 126 ierr = PetscOptionsReal("-ts_rtol","Relative tolerance for local truncation error","TSSetTolerances",ts->rtol,&ts->rtol,NULL);CHKERRQ(ierr); 127 ierr = PetscOptionsReal("-ts_atol","Absolute tolerance for local truncation error","TSSetTolerances",ts->atol,&ts->atol,NULL);CHKERRQ(ierr); 128 129 ierr = PetscOptionsBool("-ts_rhs_jacobian_test_mult","Test the RHS Jacobian for consistency with RHS at each solve ","None",ts->testjacobian,&ts->testjacobian,NULL);CHKERRQ(ierr); 130 ierr = PetscOptionsBool("-ts_rhs_jacobian_test_mult_transpose","Test the RHS Jacobian transpose for consistency with RHS at each solve ","None",ts->testjacobiantranspose,&ts->testjacobiantranspose,NULL);CHKERRQ(ierr); 131 ierr = PetscOptionsBool("-ts_use_splitrhsfunction","Use the split RHS function for multirate solvers ","TSSetUseSplitRHSFunction",ts->use_splitrhsfunction,&ts->use_splitrhsfunction,NULL);CHKERRQ(ierr); 132 #if defined(PETSC_HAVE_SAWS) 133 { 134 PetscBool set; 135 flg = PETSC_FALSE; 136 ierr = PetscOptionsBool("-ts_saws_block","Block for SAWs memory snooper at end of TSSolve","PetscObjectSAWsBlock",((PetscObject)ts)->amspublishblock,&flg,&set);CHKERRQ(ierr); 137 if (set) { 138 ierr = PetscObjectSAWsSetBlock((PetscObject)ts,flg);CHKERRQ(ierr); 139 } 140 } 141 #endif 142 143 /* Monitor options */ 144 ierr = PetscOptionsInt("-ts_monitor_frequency", "Number of time steps between monitor output", "TSMonitorSetFrequency", ts->monitorFrequency, &ts->monitorFrequency, NULL);CHKERRQ(ierr); 145 ierr = TSMonitorSetFromOptions(ts,"-ts_monitor","Monitor time and timestep size","TSMonitorDefault",TSMonitorDefault,NULL);CHKERRQ(ierr); 146 ierr = TSMonitorSetFromOptions(ts,"-ts_monitor_extreme","Monitor extreme values of the solution","TSMonitorExtreme",TSMonitorExtreme,NULL);CHKERRQ(ierr); 147 ierr = TSMonitorSetFromOptions(ts,"-ts_monitor_solution","View the solution at each timestep","TSMonitorSolution",TSMonitorSolution,NULL);CHKERRQ(ierr); 148 ierr = TSMonitorSetFromOptions(ts,"-ts_dmswarm_monitor_moments","Monitor moments of particle distribution","TSDMSwarmMonitorMoments",TSDMSwarmMonitorMoments,NULL);CHKERRQ(ierr); 149 150 ierr = PetscOptionsString("-ts_monitor_python","Use Python function","TSMonitorSet",NULL,monfilename,sizeof(monfilename),&flg);CHKERRQ(ierr); 151 if (flg) {ierr = PetscPythonMonitorSet((PetscObject)ts,monfilename);CHKERRQ(ierr);} 152 153 ierr = PetscOptionsName("-ts_monitor_lg_solution","Monitor solution graphically","TSMonitorLGSolution",&opt);CHKERRQ(ierr); 154 if (opt) { 155 PetscInt howoften = 1; 156 DM dm; 157 PetscBool net; 158 159 ierr = PetscOptionsInt("-ts_monitor_lg_solution","Monitor solution graphically","TSMonitorLGSolution",howoften,&howoften,NULL);CHKERRQ(ierr); 160 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 161 ierr = PetscObjectTypeCompare((PetscObject)dm,DMNETWORK,&net);CHKERRQ(ierr); 162 if (net) { 163 TSMonitorLGCtxNetwork ctx; 164 ierr = TSMonitorLGCtxNetworkCreate(ts,NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&ctx);CHKERRQ(ierr); 165 ierr = TSMonitorSet(ts,TSMonitorLGCtxNetworkSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxNetworkDestroy);CHKERRQ(ierr); 166 ierr = PetscOptionsBool("-ts_monitor_lg_solution_semilogy","Plot the solution with a semi-log axis","",ctx->semilogy,&ctx->semilogy,NULL);CHKERRQ(ierr); 167 } else { 168 TSMonitorLGCtx ctx; 169 ierr = TSMonitorLGCtxCreate(PETSC_COMM_SELF,NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);CHKERRQ(ierr); 170 ierr = TSMonitorSet(ts,TSMonitorLGSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);CHKERRQ(ierr); 171 } 172 } 173 174 ierr = PetscOptionsName("-ts_monitor_lg_error","Monitor error graphically","TSMonitorLGError",&opt);CHKERRQ(ierr); 175 if (opt) { 176 TSMonitorLGCtx ctx; 177 PetscInt howoften = 1; 178 179 ierr = PetscOptionsInt("-ts_monitor_lg_error","Monitor error graphically","TSMonitorLGError",howoften,&howoften,NULL);CHKERRQ(ierr); 180 ierr = TSMonitorLGCtxCreate(PETSC_COMM_SELF,NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);CHKERRQ(ierr); 181 ierr = TSMonitorSet(ts,TSMonitorLGError,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);CHKERRQ(ierr); 182 } 183 ierr = TSMonitorSetFromOptions(ts,"-ts_monitor_error","View the error at each timestep","TSMonitorError",TSMonitorError,NULL);CHKERRQ(ierr); 184 185 ierr = PetscOptionsName("-ts_monitor_lg_timestep","Monitor timestep size graphically","TSMonitorLGTimeStep",&opt);CHKERRQ(ierr); 186 if (opt) { 187 TSMonitorLGCtx ctx; 188 PetscInt howoften = 1; 189 190 ierr = PetscOptionsInt("-ts_monitor_lg_timestep","Monitor timestep size graphically","TSMonitorLGTimeStep",howoften,&howoften,NULL);CHKERRQ(ierr); 191 ierr = TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);CHKERRQ(ierr); 192 ierr = TSMonitorSet(ts,TSMonitorLGTimeStep,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);CHKERRQ(ierr); 193 } 194 ierr = PetscOptionsName("-ts_monitor_lg_timestep_log","Monitor log timestep size graphically","TSMonitorLGTimeStep",&opt);CHKERRQ(ierr); 195 if (opt) { 196 TSMonitorLGCtx ctx; 197 PetscInt howoften = 1; 198 199 ierr = PetscOptionsInt("-ts_monitor_lg_timestep_log","Monitor log timestep size graphically","TSMonitorLGTimeStep",howoften,&howoften,NULL);CHKERRQ(ierr); 200 ierr = TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);CHKERRQ(ierr); 201 ierr = TSMonitorSet(ts,TSMonitorLGTimeStep,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);CHKERRQ(ierr); 202 ctx->semilogy = PETSC_TRUE; 203 } 204 205 ierr = PetscOptionsName("-ts_monitor_lg_snes_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGSNESIterations",&opt);CHKERRQ(ierr); 206 if (opt) { 207 TSMonitorLGCtx ctx; 208 PetscInt howoften = 1; 209 210 ierr = PetscOptionsInt("-ts_monitor_lg_snes_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGSNESIterations",howoften,&howoften,NULL);CHKERRQ(ierr); 211 ierr = TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);CHKERRQ(ierr); 212 ierr = TSMonitorSet(ts,TSMonitorLGSNESIterations,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);CHKERRQ(ierr); 213 } 214 ierr = PetscOptionsName("-ts_monitor_lg_ksp_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGKSPIterations",&opt);CHKERRQ(ierr); 215 if (opt) { 216 TSMonitorLGCtx ctx; 217 PetscInt howoften = 1; 218 219 ierr = PetscOptionsInt("-ts_monitor_lg_ksp_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGKSPIterations",howoften,&howoften,NULL);CHKERRQ(ierr); 220 ierr = TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);CHKERRQ(ierr); 221 ierr = TSMonitorSet(ts,TSMonitorLGKSPIterations,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);CHKERRQ(ierr); 222 } 223 ierr = PetscOptionsName("-ts_monitor_sp_eig","Monitor eigenvalues of linearized operator graphically","TSMonitorSPEig",&opt);CHKERRQ(ierr); 224 if (opt) { 225 TSMonitorSPEigCtx ctx; 226 PetscInt howoften = 1; 227 228 ierr = PetscOptionsInt("-ts_monitor_sp_eig","Monitor eigenvalues of linearized operator graphically","TSMonitorSPEig",howoften,&howoften,NULL);CHKERRQ(ierr); 229 ierr = TSMonitorSPEigCtxCreate(PETSC_COMM_SELF,NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);CHKERRQ(ierr); 230 ierr = TSMonitorSet(ts,TSMonitorSPEig,ctx,(PetscErrorCode (*)(void**))TSMonitorSPEigCtxDestroy);CHKERRQ(ierr); 231 } 232 ierr = PetscOptionsName("-ts_monitor_sp_swarm","Display particle phase from the DMSwarm","TSMonitorSPSwarm",&opt);CHKERRQ(ierr); 233 if (opt) { 234 TSMonitorSPCtx ctx; 235 PetscInt howoften = 1; 236 ierr = PetscOptionsInt("-ts_monitor_sp_swarm","Display particles phase from the DMSwarm","TSMonitorSPSwarm",howoften,&howoften,NULL);CHKERRQ(ierr); 237 ierr = TSMonitorSPCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx);CHKERRQ(ierr); 238 ierr = TSMonitorSet(ts, TSMonitorSPSwarmSolution, ctx, (PetscErrorCode (*)(void**))TSMonitorSPCtxDestroy);CHKERRQ(ierr); 239 } 240 opt = PETSC_FALSE; 241 ierr = PetscOptionsName("-ts_monitor_draw_solution","Monitor solution graphically","TSMonitorDrawSolution",&opt);CHKERRQ(ierr); 242 if (opt) { 243 TSMonitorDrawCtx ctx; 244 PetscInt howoften = 1; 245 246 ierr = PetscOptionsInt("-ts_monitor_draw_solution","Monitor solution graphically","TSMonitorDrawSolution",howoften,&howoften,NULL);CHKERRQ(ierr); 247 ierr = TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),NULL,"Computed Solution",PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);CHKERRQ(ierr); 248 ierr = TSMonitorSet(ts,TSMonitorDrawSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);CHKERRQ(ierr); 249 } 250 opt = PETSC_FALSE; 251 ierr = PetscOptionsName("-ts_monitor_draw_solution_phase","Monitor solution graphically","TSMonitorDrawSolutionPhase",&opt);CHKERRQ(ierr); 252 if (opt) { 253 TSMonitorDrawCtx ctx; 254 PetscReal bounds[4]; 255 PetscInt n = 4; 256 PetscDraw draw; 257 PetscDrawAxis axis; 258 259 ierr = PetscOptionsRealArray("-ts_monitor_draw_solution_phase","Monitor solution graphically","TSMonitorDrawSolutionPhase",bounds,&n,NULL);CHKERRQ(ierr); 260 if (n != 4) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Must provide bounding box of phase field"); 261 ierr = TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,300,300,1,&ctx);CHKERRQ(ierr); 262 ierr = PetscViewerDrawGetDraw(ctx->viewer,0,&draw);CHKERRQ(ierr); 263 ierr = PetscViewerDrawGetDrawAxis(ctx->viewer,0,&axis);CHKERRQ(ierr); 264 ierr = PetscDrawAxisSetLimits(axis,bounds[0],bounds[2],bounds[1],bounds[3]);CHKERRQ(ierr); 265 ierr = PetscDrawAxisSetLabels(axis,"Phase Diagram","Variable 1","Variable 2");CHKERRQ(ierr); 266 ierr = TSMonitorSet(ts,TSMonitorDrawSolutionPhase,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);CHKERRQ(ierr); 267 } 268 opt = PETSC_FALSE; 269 ierr = PetscOptionsName("-ts_monitor_draw_error","Monitor error graphically","TSMonitorDrawError",&opt);CHKERRQ(ierr); 270 if (opt) { 271 TSMonitorDrawCtx ctx; 272 PetscInt howoften = 1; 273 274 ierr = PetscOptionsInt("-ts_monitor_draw_error","Monitor error graphically","TSMonitorDrawError",howoften,&howoften,NULL);CHKERRQ(ierr); 275 ierr = TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),NULL,"Error",PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);CHKERRQ(ierr); 276 ierr = TSMonitorSet(ts,TSMonitorDrawError,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);CHKERRQ(ierr); 277 } 278 opt = PETSC_FALSE; 279 ierr = PetscOptionsName("-ts_monitor_draw_solution_function","Monitor solution provided by TSMonitorSetSolutionFunction() graphically","TSMonitorDrawSolutionFunction",&opt);CHKERRQ(ierr); 280 if (opt) { 281 TSMonitorDrawCtx ctx; 282 PetscInt howoften = 1; 283 284 ierr = PetscOptionsInt("-ts_monitor_draw_solution_function","Monitor solution provided by TSMonitorSetSolutionFunction() graphically","TSMonitorDrawSolutionFunction",howoften,&howoften,NULL);CHKERRQ(ierr); 285 ierr = TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),NULL,"Solution provided by user function",PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);CHKERRQ(ierr); 286 ierr = TSMonitorSet(ts,TSMonitorDrawSolutionFunction,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);CHKERRQ(ierr); 287 } 288 289 opt = PETSC_FALSE; 290 ierr = PetscOptionsString("-ts_monitor_solution_vtk","Save each time step to a binary file, use filename-%%03D.vts","TSMonitorSolutionVTK",NULL,monfilename,sizeof(monfilename),&flg);CHKERRQ(ierr); 291 if (flg) { 292 const char *ptr,*ptr2; 293 char *filetemplate; 294 if (!monfilename[0]) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts"); 295 /* Do some cursory validation of the input. */ 296 ierr = PetscStrstr(monfilename,"%",(char**)&ptr);CHKERRQ(ierr); 297 if (!ptr) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts"); 298 for (ptr++; ptr && *ptr; ptr++) { 299 ierr = PetscStrchr("DdiouxX",*ptr,(char**)&ptr2);CHKERRQ(ierr); 300 if (!ptr2 && (*ptr < '0' || '9' < *ptr)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Invalid file template argument to -ts_monitor_solution_vtk, should look like filename-%%03D.vts"); 301 if (ptr2) break; 302 } 303 ierr = PetscStrallocpy(monfilename,&filetemplate);CHKERRQ(ierr); 304 ierr = TSMonitorSet(ts,TSMonitorSolutionVTK,filetemplate,(PetscErrorCode (*)(void**))TSMonitorSolutionVTKDestroy);CHKERRQ(ierr); 305 } 306 307 ierr = PetscOptionsString("-ts_monitor_dmda_ray","Display a ray of the solution","None","y=0",dir,sizeof(dir),&flg);CHKERRQ(ierr); 308 if (flg) { 309 TSMonitorDMDARayCtx *rayctx; 310 int ray = 0; 311 DMDirection ddir; 312 DM da; 313 PetscMPIInt rank; 314 315 if (dir[1] != '=') SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Unknown ray %s",dir); 316 if (dir[0] == 'x') ddir = DM_X; 317 else if (dir[0] == 'y') ddir = DM_Y; 318 else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Unknown ray %s",dir); 319 sscanf(dir+2,"%d",&ray); 320 321 ierr = PetscInfo2(((PetscObject)ts),"Displaying DMDA ray %c = %d\n",dir[0],ray);CHKERRQ(ierr); 322 ierr = PetscNew(&rayctx);CHKERRQ(ierr); 323 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 324 ierr = DMDAGetRay(da,ddir,ray,&rayctx->ray,&rayctx->scatter);CHKERRQ(ierr); 325 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)ts),&rank);CHKERRMPI(ierr); 326 if (!rank) { 327 ierr = PetscViewerDrawOpen(PETSC_COMM_SELF,NULL,NULL,0,0,600,300,&rayctx->viewer);CHKERRQ(ierr); 328 } 329 rayctx->lgctx = NULL; 330 ierr = TSMonitorSet(ts,TSMonitorDMDARay,rayctx,TSMonitorDMDARayDestroy);CHKERRQ(ierr); 331 } 332 ierr = PetscOptionsString("-ts_monitor_lg_dmda_ray","Display a ray of the solution","None","x=0",dir,sizeof(dir),&flg);CHKERRQ(ierr); 333 if (flg) { 334 TSMonitorDMDARayCtx *rayctx; 335 int ray = 0; 336 DMDirection ddir; 337 DM da; 338 PetscInt howoften = 1; 339 340 if (dir[1] != '=') SETERRQ1(PetscObjectComm((PetscObject) ts), PETSC_ERR_ARG_WRONG, "Malformed ray %s", dir); 341 if (dir[0] == 'x') ddir = DM_X; 342 else if (dir[0] == 'y') ddir = DM_Y; 343 else SETERRQ1(PetscObjectComm((PetscObject) ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir); 344 sscanf(dir+2, "%d", &ray); 345 346 ierr = PetscInfo2(((PetscObject) ts),"Displaying LG DMDA ray %c = %d\n", dir[0], ray);CHKERRQ(ierr); 347 ierr = PetscNew(&rayctx);CHKERRQ(ierr); 348 ierr = TSGetDM(ts, &da);CHKERRQ(ierr); 349 ierr = DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter);CHKERRQ(ierr); 350 ierr = TSMonitorLGCtxCreate(PETSC_COMM_SELF,NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&rayctx->lgctx);CHKERRQ(ierr); 351 ierr = TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy);CHKERRQ(ierr); 352 } 353 354 ierr = PetscOptionsName("-ts_monitor_envelope","Monitor maximum and minimum value of each component of the solution","TSMonitorEnvelope",&opt);CHKERRQ(ierr); 355 if (opt) { 356 TSMonitorEnvelopeCtx ctx; 357 358 ierr = TSMonitorEnvelopeCtxCreate(ts,&ctx);CHKERRQ(ierr); 359 ierr = TSMonitorSet(ts,TSMonitorEnvelope,ctx,(PetscErrorCode (*)(void**))TSMonitorEnvelopeCtxDestroy);CHKERRQ(ierr); 360 } 361 flg = PETSC_FALSE; 362 ierr = PetscOptionsBool("-ts_monitor_cancel","Remove all monitors","TSMonitorCancel",flg,&flg,&opt);CHKERRQ(ierr); 363 if (opt && flg) {ierr = TSMonitorCancel(ts);CHKERRQ(ierr);} 364 365 flg = PETSC_FALSE; 366 ierr = PetscOptionsBool("-ts_fd_color", "Use finite differences with coloring to compute IJacobian", "TSComputeJacobianDefaultColor", flg, &flg, NULL);CHKERRQ(ierr); 367 if (flg) { 368 DM dm; 369 DMTS tdm; 370 371 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 372 ierr = DMGetDMTS(dm, &tdm);CHKERRQ(ierr); 373 tdm->ijacobianctx = NULL; 374 ierr = TSSetIJacobian(ts, NULL, NULL, TSComputeIJacobianDefaultColor, NULL);CHKERRQ(ierr); 375 ierr = PetscInfo(ts, "Setting default finite difference coloring Jacobian matrix\n");CHKERRQ(ierr); 376 } 377 378 /* Handle specific TS options */ 379 if (ts->ops->setfromoptions) { 380 ierr = (*ts->ops->setfromoptions)(PetscOptionsObject,ts);CHKERRQ(ierr); 381 } 382 383 /* Handle TSAdapt options */ 384 ierr = TSGetAdapt(ts,&ts->adapt);CHKERRQ(ierr); 385 ierr = TSAdaptSetDefaultType(ts->adapt,ts->default_adapt_type);CHKERRQ(ierr); 386 ierr = TSAdaptSetFromOptions(PetscOptionsObject,ts->adapt);CHKERRQ(ierr); 387 388 /* TS trajectory must be set after TS, since it may use some TS options above */ 389 tflg = ts->trajectory ? PETSC_TRUE : PETSC_FALSE; 390 ierr = PetscOptionsBool("-ts_save_trajectory","Save the solution at each timestep","TSSetSaveTrajectory",tflg,&tflg,NULL);CHKERRQ(ierr); 391 if (tflg) { 392 ierr = TSSetSaveTrajectory(ts);CHKERRQ(ierr); 393 } 394 395 ierr = TSAdjointSetFromOptions(PetscOptionsObject,ts);CHKERRQ(ierr); 396 397 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 398 ierr = PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)ts);CHKERRQ(ierr); 399 ierr = PetscOptionsEnd();CHKERRQ(ierr); 400 401 if (ts->trajectory) { 402 ierr = TSTrajectorySetFromOptions(ts->trajectory,ts);CHKERRQ(ierr); 403 } 404 405 /* why do we have to do this here and not during TSSetUp? */ 406 ierr = TSGetSNES(ts,&ts->snes);CHKERRQ(ierr); 407 if (ts->problem_type == TS_LINEAR) { 408 ierr = PetscObjectTypeCompareAny((PetscObject)ts->snes,&flg,SNESKSPONLY,SNESKSPTRANSPOSEONLY,"");CHKERRQ(ierr); 409 if (!flg) { ierr = SNESSetType(ts->snes,SNESKSPONLY);CHKERRQ(ierr); } 410 } 411 ierr = SNESSetFromOptions(ts->snes);CHKERRQ(ierr); 412 PetscFunctionReturn(0); 413 } 414 415 /*@ 416 TSGetTrajectory - Gets the trajectory from a TS if it exists 417 418 Collective on TS 419 420 Input Parameters: 421 . ts - the TS context obtained from TSCreate() 422 423 Output Parameters: 424 . tr - the TSTrajectory object, if it exists 425 426 Note: This routine should be called after all TS options have been set 427 428 Level: advanced 429 430 .seealso: TSGetTrajectory(), TSAdjointSolve(), TSTrajectory, TSTrajectoryCreate() 431 432 @*/ 433 PetscErrorCode TSGetTrajectory(TS ts,TSTrajectory *tr) 434 { 435 PetscFunctionBegin; 436 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 437 *tr = ts->trajectory; 438 PetscFunctionReturn(0); 439 } 440 441 /*@ 442 TSSetSaveTrajectory - Causes the TS to save its solutions as it iterates forward in time in a TSTrajectory object 443 444 Collective on TS 445 446 Input Parameters: 447 . ts - the TS context obtained from TSCreate() 448 449 Options Database: 450 + -ts_save_trajectory - saves the trajectory to a file 451 - -ts_trajectory_type type 452 453 Note: This routine should be called after all TS options have been set 454 455 The TSTRAJECTORYVISUALIZATION files can be loaded into Python with $PETSC_DIR/lib/petsc/bin/PetscBinaryIOTrajectory.py and 456 MATLAB with $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m 457 458 Level: intermediate 459 460 .seealso: TSGetTrajectory(), TSAdjointSolve() 461 462 @*/ 463 PetscErrorCode TSSetSaveTrajectory(TS ts) 464 { 465 PetscErrorCode ierr; 466 467 PetscFunctionBegin; 468 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 469 if (!ts->trajectory) { 470 ierr = TSTrajectoryCreate(PetscObjectComm((PetscObject)ts),&ts->trajectory);CHKERRQ(ierr); 471 } 472 PetscFunctionReturn(0); 473 } 474 475 /*@ 476 TSResetTrajectory - Destroys and recreates the internal TSTrajectory object 477 478 Collective on TS 479 480 Input Parameters: 481 . ts - the TS context obtained from TSCreate() 482 483 Level: intermediate 484 485 .seealso: TSGetTrajectory(), TSAdjointSolve() 486 487 @*/ 488 PetscErrorCode TSResetTrajectory(TS ts) 489 { 490 PetscErrorCode ierr; 491 492 PetscFunctionBegin; 493 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 494 if (ts->trajectory) { 495 ierr = TSTrajectoryDestroy(&ts->trajectory);CHKERRQ(ierr); 496 ierr = TSTrajectoryCreate(PetscObjectComm((PetscObject)ts),&ts->trajectory);CHKERRQ(ierr); 497 } 498 PetscFunctionReturn(0); 499 } 500 501 /*@ 502 TSComputeRHSJacobian - Computes the Jacobian matrix that has been 503 set with TSSetRHSJacobian(). 504 505 Collective on TS 506 507 Input Parameters: 508 + ts - the TS context 509 . t - current timestep 510 - U - input vector 511 512 Output Parameters: 513 + A - Jacobian matrix 514 . B - optional preconditioning matrix 515 - flag - flag indicating matrix structure 516 517 Notes: 518 Most users should not need to explicitly call this routine, as it 519 is used internally within the nonlinear solvers. 520 521 See KSPSetOperators() for important information about setting the 522 flag parameter. 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/contex to compute it. 1291 1292 Not Collective 1293 1294 Input Parameter: 1295 . ts - the TS context 1296 1297 Output Parameter: 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 Parameter: 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 Arguments: 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/contex to compute it. 1485 1486 Not Collective 1487 1488 Input Parameter: 1489 . ts - the TS context 1490 1491 Output Parameter: 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 Arguments: 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 Parameter: 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) { 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) { 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 Parameter: 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->num_rhs_splits = 0; 2804 ts->setupcalled = PETSC_FALSE; 2805 PetscFunctionReturn(0); 2806 } 2807 2808 /*@C 2809 TSDestroy - Destroys the timestepper context that was created 2810 with TSCreate(). 2811 2812 Collective on TS 2813 2814 Input Parameter: 2815 . ts - the TS context obtained from TSCreate() 2816 2817 Level: beginner 2818 2819 .seealso: TSCreate(), TSSetUp(), TSSolve() 2820 @*/ 2821 PetscErrorCode TSDestroy(TS *ts) 2822 { 2823 PetscErrorCode ierr; 2824 2825 PetscFunctionBegin; 2826 if (!*ts) PetscFunctionReturn(0); 2827 PetscValidHeaderSpecific(*ts,TS_CLASSID,1); 2828 if (--((PetscObject)(*ts))->refct > 0) {*ts = NULL; PetscFunctionReturn(0);} 2829 2830 ierr = TSReset(*ts);CHKERRQ(ierr); 2831 ierr = TSAdjointReset(*ts);CHKERRQ(ierr); 2832 if ((*ts)->forward_solve) { 2833 ierr = TSForwardReset(*ts);CHKERRQ(ierr); 2834 } 2835 /* if memory was published with SAWs then destroy it */ 2836 ierr = PetscObjectSAWsViewOff((PetscObject)*ts);CHKERRQ(ierr); 2837 if ((*ts)->ops->destroy) {ierr = (*(*ts)->ops->destroy)((*ts));CHKERRQ(ierr);} 2838 2839 ierr = TSTrajectoryDestroy(&(*ts)->trajectory);CHKERRQ(ierr); 2840 2841 ierr = TSAdaptDestroy(&(*ts)->adapt);CHKERRQ(ierr); 2842 ierr = TSEventDestroy(&(*ts)->event);CHKERRQ(ierr); 2843 2844 ierr = SNESDestroy(&(*ts)->snes);CHKERRQ(ierr); 2845 ierr = DMDestroy(&(*ts)->dm);CHKERRQ(ierr); 2846 ierr = TSMonitorCancel((*ts));CHKERRQ(ierr); 2847 ierr = TSAdjointMonitorCancel((*ts));CHKERRQ(ierr); 2848 2849 ierr = TSDestroy(&(*ts)->quadraturets);CHKERRQ(ierr); 2850 ierr = PetscHeaderDestroy(ts);CHKERRQ(ierr); 2851 PetscFunctionReturn(0); 2852 } 2853 2854 /*@ 2855 TSGetSNES - Returns the SNES (nonlinear solver) associated with 2856 a TS (timestepper) context. Valid only for nonlinear problems. 2857 2858 Not Collective, but SNES is parallel if TS is parallel 2859 2860 Input Parameter: 2861 . ts - the TS context obtained from TSCreate() 2862 2863 Output Parameter: 2864 . snes - the nonlinear solver context 2865 2866 Notes: 2867 The user can then directly manipulate the SNES context to set various 2868 options, etc. Likewise, the user can then extract and manipulate the 2869 KSP, KSP, and PC contexts as well. 2870 2871 TSGetSNES() does not work for integrators that do not use SNES; in 2872 this case TSGetSNES() returns NULL in snes. 2873 2874 Level: beginner 2875 2876 @*/ 2877 PetscErrorCode TSGetSNES(TS ts,SNES *snes) 2878 { 2879 PetscErrorCode ierr; 2880 2881 PetscFunctionBegin; 2882 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2883 PetscValidPointer(snes,2); 2884 if (!ts->snes) { 2885 ierr = SNESCreate(PetscObjectComm((PetscObject)ts),&ts->snes);CHKERRQ(ierr); 2886 ierr = PetscObjectSetOptions((PetscObject)ts->snes,((PetscObject)ts)->options);CHKERRQ(ierr); 2887 ierr = SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 2888 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->snes);CHKERRQ(ierr); 2889 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);CHKERRQ(ierr); 2890 if (ts->dm) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 2891 if (ts->problem_type == TS_LINEAR) { 2892 ierr = SNESSetType(ts->snes,SNESKSPONLY);CHKERRQ(ierr); 2893 } 2894 } 2895 *snes = ts->snes; 2896 PetscFunctionReturn(0); 2897 } 2898 2899 /*@ 2900 TSSetSNES - Set the SNES (nonlinear solver) to be used by the timestepping context 2901 2902 Collective 2903 2904 Input Parameter: 2905 + ts - the TS context obtained from TSCreate() 2906 - snes - the nonlinear solver context 2907 2908 Notes: 2909 Most users should have the TS created by calling TSGetSNES() 2910 2911 Level: developer 2912 2913 @*/ 2914 PetscErrorCode TSSetSNES(TS ts,SNES snes) 2915 { 2916 PetscErrorCode ierr; 2917 PetscErrorCode (*func)(SNES,Vec,Mat,Mat,void*); 2918 2919 PetscFunctionBegin; 2920 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2921 PetscValidHeaderSpecific(snes,SNES_CLASSID,2); 2922 ierr = PetscObjectReference((PetscObject)snes);CHKERRQ(ierr); 2923 ierr = SNESDestroy(&ts->snes);CHKERRQ(ierr); 2924 2925 ts->snes = snes; 2926 2927 ierr = SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 2928 ierr = SNESGetJacobian(ts->snes,NULL,NULL,&func,NULL);CHKERRQ(ierr); 2929 if (func == SNESTSFormJacobian) { 2930 ierr = SNESSetJacobian(ts->snes,NULL,NULL,SNESTSFormJacobian,ts);CHKERRQ(ierr); 2931 } 2932 PetscFunctionReturn(0); 2933 } 2934 2935 /*@ 2936 TSGetKSP - Returns the KSP (linear solver) associated with 2937 a TS (timestepper) context. 2938 2939 Not Collective, but KSP is parallel if TS is parallel 2940 2941 Input Parameter: 2942 . ts - the TS context obtained from TSCreate() 2943 2944 Output Parameter: 2945 . ksp - the nonlinear solver context 2946 2947 Notes: 2948 The user can then directly manipulate the KSP context to set various 2949 options, etc. Likewise, the user can then extract and manipulate the 2950 KSP and PC contexts as well. 2951 2952 TSGetKSP() does not work for integrators that do not use KSP; 2953 in this case TSGetKSP() returns NULL in ksp. 2954 2955 Level: beginner 2956 2957 @*/ 2958 PetscErrorCode TSGetKSP(TS ts,KSP *ksp) 2959 { 2960 PetscErrorCode ierr; 2961 SNES snes; 2962 2963 PetscFunctionBegin; 2964 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2965 PetscValidPointer(ksp,2); 2966 if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first"); 2967 if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()"); 2968 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2969 ierr = SNESGetKSP(snes,ksp);CHKERRQ(ierr); 2970 PetscFunctionReturn(0); 2971 } 2972 2973 /* ----------- Routines to set solver parameters ---------- */ 2974 2975 /*@ 2976 TSSetMaxSteps - Sets the maximum number of steps to use. 2977 2978 Logically Collective on TS 2979 2980 Input Parameters: 2981 + ts - the TS context obtained from TSCreate() 2982 - maxsteps - maximum number of steps to use 2983 2984 Options Database Keys: 2985 . -ts_max_steps <maxsteps> - Sets maxsteps 2986 2987 Notes: 2988 The default maximum number of steps is 5000 2989 2990 Level: intermediate 2991 2992 .seealso: TSGetMaxSteps(), TSSetMaxTime(), TSSetExactFinalTime() 2993 @*/ 2994 PetscErrorCode TSSetMaxSteps(TS ts,PetscInt maxsteps) 2995 { 2996 PetscFunctionBegin; 2997 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2998 PetscValidLogicalCollectiveInt(ts,maxsteps,2); 2999 if (maxsteps < 0) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of steps must be non-negative"); 3000 ts->max_steps = maxsteps; 3001 PetscFunctionReturn(0); 3002 } 3003 3004 /*@ 3005 TSGetMaxSteps - Gets the maximum number of steps to use. 3006 3007 Not Collective 3008 3009 Input Parameters: 3010 . ts - the TS context obtained from TSCreate() 3011 3012 Output Parameter: 3013 . maxsteps - maximum number of steps to use 3014 3015 Level: advanced 3016 3017 .seealso: TSSetMaxSteps(), TSGetMaxTime(), TSSetMaxTime() 3018 @*/ 3019 PetscErrorCode TSGetMaxSteps(TS ts,PetscInt *maxsteps) 3020 { 3021 PetscFunctionBegin; 3022 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3023 PetscValidIntPointer(maxsteps,2); 3024 *maxsteps = ts->max_steps; 3025 PetscFunctionReturn(0); 3026 } 3027 3028 /*@ 3029 TSSetMaxTime - Sets the maximum (or final) time for timestepping. 3030 3031 Logically Collective on TS 3032 3033 Input Parameters: 3034 + ts - the TS context obtained from TSCreate() 3035 - maxtime - final time to step to 3036 3037 Options Database Keys: 3038 . -ts_max_time <maxtime> - Sets maxtime 3039 3040 Notes: 3041 The default maximum time is 5.0 3042 3043 Level: intermediate 3044 3045 .seealso: TSGetMaxTime(), TSSetMaxSteps(), TSSetExactFinalTime() 3046 @*/ 3047 PetscErrorCode TSSetMaxTime(TS ts,PetscReal maxtime) 3048 { 3049 PetscFunctionBegin; 3050 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3051 PetscValidLogicalCollectiveReal(ts,maxtime,2); 3052 ts->max_time = maxtime; 3053 PetscFunctionReturn(0); 3054 } 3055 3056 /*@ 3057 TSGetMaxTime - Gets the maximum (or final) time for timestepping. 3058 3059 Not Collective 3060 3061 Input Parameters: 3062 . ts - the TS context obtained from TSCreate() 3063 3064 Output Parameter: 3065 . maxtime - final time to step to 3066 3067 Level: advanced 3068 3069 .seealso: TSSetMaxTime(), TSGetMaxSteps(), TSSetMaxSteps() 3070 @*/ 3071 PetscErrorCode TSGetMaxTime(TS ts,PetscReal *maxtime) 3072 { 3073 PetscFunctionBegin; 3074 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3075 PetscValidRealPointer(maxtime,2); 3076 *maxtime = ts->max_time; 3077 PetscFunctionReturn(0); 3078 } 3079 3080 /*@ 3081 TSSetInitialTimeStep - Deprecated, use TSSetTime() and TSSetTimeStep(). 3082 3083 Level: deprecated 3084 3085 @*/ 3086 PetscErrorCode TSSetInitialTimeStep(TS ts,PetscReal initial_time,PetscReal time_step) 3087 { 3088 PetscErrorCode ierr; 3089 PetscFunctionBegin; 3090 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3091 ierr = TSSetTime(ts,initial_time);CHKERRQ(ierr); 3092 ierr = TSSetTimeStep(ts,time_step);CHKERRQ(ierr); 3093 PetscFunctionReturn(0); 3094 } 3095 3096 /*@ 3097 TSGetDuration - Deprecated, use TSGetMaxSteps() and TSGetMaxTime(). 3098 3099 Level: deprecated 3100 3101 @*/ 3102 PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime) 3103 { 3104 PetscFunctionBegin; 3105 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 3106 if (maxsteps) { 3107 PetscValidIntPointer(maxsteps,2); 3108 *maxsteps = ts->max_steps; 3109 } 3110 if (maxtime) { 3111 PetscValidRealPointer(maxtime,3); 3112 *maxtime = ts->max_time; 3113 } 3114 PetscFunctionReturn(0); 3115 } 3116 3117 /*@ 3118 TSSetDuration - Deprecated, use TSSetMaxSteps() and TSSetMaxTime(). 3119 3120 Level: deprecated 3121 3122 @*/ 3123 PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime) 3124 { 3125 PetscFunctionBegin; 3126 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3127 PetscValidLogicalCollectiveInt(ts,maxsteps,2); 3128 PetscValidLogicalCollectiveReal(ts,maxtime,3); 3129 if (maxsteps >= 0) ts->max_steps = maxsteps; 3130 if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime; 3131 PetscFunctionReturn(0); 3132 } 3133 3134 /*@ 3135 TSGetTimeStepNumber - Deprecated, use TSGetStepNumber(). 3136 3137 Level: deprecated 3138 3139 @*/ 3140 PetscErrorCode TSGetTimeStepNumber(TS ts,PetscInt *steps) { return TSGetStepNumber(ts,steps); } 3141 3142 /*@ 3143 TSGetTotalSteps - Deprecated, use TSGetStepNumber(). 3144 3145 Level: deprecated 3146 3147 @*/ 3148 PetscErrorCode TSGetTotalSteps(TS ts,PetscInt *steps) { return TSGetStepNumber(ts,steps); } 3149 3150 /*@ 3151 TSSetSolution - Sets the initial solution vector 3152 for use by the TS routines. 3153 3154 Logically Collective on TS 3155 3156 Input Parameters: 3157 + ts - the TS context obtained from TSCreate() 3158 - u - the solution vector 3159 3160 Level: beginner 3161 3162 .seealso: TSSetSolutionFunction(), TSGetSolution(), TSCreate() 3163 @*/ 3164 PetscErrorCode TSSetSolution(TS ts,Vec u) 3165 { 3166 PetscErrorCode ierr; 3167 DM dm; 3168 3169 PetscFunctionBegin; 3170 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3171 PetscValidHeaderSpecific(u,VEC_CLASSID,2); 3172 ierr = PetscObjectReference((PetscObject)u);CHKERRQ(ierr); 3173 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 3174 ts->vec_sol = u; 3175 3176 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 3177 ierr = DMShellSetGlobalVector(dm,u);CHKERRQ(ierr); 3178 PetscFunctionReturn(0); 3179 } 3180 3181 /*@C 3182 TSSetPreStep - Sets the general-purpose function 3183 called once at the beginning of each time step. 3184 3185 Logically Collective on TS 3186 3187 Input Parameters: 3188 + ts - The TS context obtained from TSCreate() 3189 - func - The function 3190 3191 Calling sequence of func: 3192 . PetscErrorCode func (TS ts); 3193 3194 Level: intermediate 3195 3196 .seealso: TSSetPreStage(), TSSetPostStage(), TSSetPostStep(), TSStep(), TSRestartStep() 3197 @*/ 3198 PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS)) 3199 { 3200 PetscFunctionBegin; 3201 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 3202 ts->prestep = func; 3203 PetscFunctionReturn(0); 3204 } 3205 3206 /*@ 3207 TSPreStep - Runs the user-defined pre-step function. 3208 3209 Collective on TS 3210 3211 Input Parameters: 3212 . ts - The TS context obtained from TSCreate() 3213 3214 Notes: 3215 TSPreStep() is typically used within time stepping implementations, 3216 so most users would not generally call this routine themselves. 3217 3218 Level: developer 3219 3220 .seealso: TSSetPreStep(), TSPreStage(), TSPostStage(), TSPostStep() 3221 @*/ 3222 PetscErrorCode TSPreStep(TS ts) 3223 { 3224 PetscErrorCode ierr; 3225 3226 PetscFunctionBegin; 3227 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3228 if (ts->prestep) { 3229 Vec U; 3230 PetscObjectState sprev,spost; 3231 3232 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 3233 ierr = PetscObjectStateGet((PetscObject)U,&sprev);CHKERRQ(ierr); 3234 PetscStackCallStandard((*ts->prestep),(ts)); 3235 ierr = PetscObjectStateGet((PetscObject)U,&spost);CHKERRQ(ierr); 3236 if (sprev != spost) {ierr = TSRestartStep(ts);CHKERRQ(ierr);} 3237 } 3238 PetscFunctionReturn(0); 3239 } 3240 3241 /*@C 3242 TSSetPreStage - Sets the general-purpose function 3243 called once at the beginning of each stage. 3244 3245 Logically Collective on TS 3246 3247 Input Parameters: 3248 + ts - The TS context obtained from TSCreate() 3249 - func - The function 3250 3251 Calling sequence of func: 3252 . PetscErrorCode func(TS ts, PetscReal stagetime); 3253 3254 Level: intermediate 3255 3256 Note: 3257 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 3258 The time step number being computed can be queried using TSGetStepNumber() and the total size of the step being 3259 attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime(). 3260 3261 .seealso: TSSetPostStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 3262 @*/ 3263 PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS,PetscReal)) 3264 { 3265 PetscFunctionBegin; 3266 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 3267 ts->prestage = func; 3268 PetscFunctionReturn(0); 3269 } 3270 3271 /*@C 3272 TSSetPostStage - Sets the general-purpose function 3273 called once at the end of each stage. 3274 3275 Logically Collective on TS 3276 3277 Input Parameters: 3278 + ts - The TS context obtained from TSCreate() 3279 - func - The function 3280 3281 Calling sequence of func: 3282 . PetscErrorCode func(TS ts, PetscReal stagetime, PetscInt stageindex, Vec* Y); 3283 3284 Level: intermediate 3285 3286 Note: 3287 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 3288 The time step number being computed can be queried using TSGetStepNumber() and the total size of the step being 3289 attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime(). 3290 3291 .seealso: TSSetPreStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 3292 @*/ 3293 PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS,PetscReal,PetscInt,Vec*)) 3294 { 3295 PetscFunctionBegin; 3296 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 3297 ts->poststage = func; 3298 PetscFunctionReturn(0); 3299 } 3300 3301 /*@C 3302 TSSetPostEvaluate - Sets the general-purpose function 3303 called once at the end of each step evaluation. 3304 3305 Logically Collective on TS 3306 3307 Input Parameters: 3308 + ts - The TS context obtained from TSCreate() 3309 - func - The function 3310 3311 Calling sequence of func: 3312 . PetscErrorCode func(TS ts); 3313 3314 Level: intermediate 3315 3316 Note: 3317 Semantically, TSSetPostEvaluate() differs from TSSetPostStep() since the function it sets is called before event-handling 3318 thus guaranteeing the same solution (computed by the time-stepper) will be passed to it. On the other hand, TSPostStep() 3319 may be passed a different solution, possibly changed by the event handler. TSPostEvaluate() is called after the next step 3320 solution is evaluated allowing to modify it, if need be. The solution can be obtained with TSGetSolution(), the time step 3321 with TSGetTimeStep(), and the time at the start of the step is available via TSGetTime() 3322 3323 .seealso: TSSetPreStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 3324 @*/ 3325 PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS)) 3326 { 3327 PetscFunctionBegin; 3328 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 3329 ts->postevaluate = func; 3330 PetscFunctionReturn(0); 3331 } 3332 3333 /*@ 3334 TSPreStage - Runs the user-defined pre-stage function set using TSSetPreStage() 3335 3336 Collective on TS 3337 3338 Input Parameters: 3339 . ts - The TS context obtained from TSCreate() 3340 stagetime - The absolute time of the current stage 3341 3342 Notes: 3343 TSPreStage() is typically used within time stepping implementations, 3344 most users would not generally call this routine themselves. 3345 3346 Level: developer 3347 3348 .seealso: TSPostStage(), TSSetPreStep(), TSPreStep(), TSPostStep() 3349 @*/ 3350 PetscErrorCode TSPreStage(TS ts, PetscReal stagetime) 3351 { 3352 PetscFunctionBegin; 3353 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3354 if (ts->prestage) { 3355 PetscStackCallStandard((*ts->prestage),(ts,stagetime)); 3356 } 3357 PetscFunctionReturn(0); 3358 } 3359 3360 /*@ 3361 TSPostStage - Runs the user-defined post-stage function set using TSSetPostStage() 3362 3363 Collective on TS 3364 3365 Input Parameters: 3366 . ts - The TS context obtained from TSCreate() 3367 stagetime - The absolute time of the current stage 3368 stageindex - Stage number 3369 Y - Array of vectors (of size = total number 3370 of stages) with the stage solutions 3371 3372 Notes: 3373 TSPostStage() is typically used within time stepping implementations, 3374 most users would not generally call this routine themselves. 3375 3376 Level: developer 3377 3378 .seealso: TSPreStage(), TSSetPreStep(), TSPreStep(), TSPostStep() 3379 @*/ 3380 PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y) 3381 { 3382 PetscFunctionBegin; 3383 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3384 if (ts->poststage) { 3385 PetscStackCallStandard((*ts->poststage),(ts,stagetime,stageindex,Y)); 3386 } 3387 PetscFunctionReturn(0); 3388 } 3389 3390 /*@ 3391 TSPostEvaluate - Runs the user-defined post-evaluate function set using TSSetPostEvaluate() 3392 3393 Collective on TS 3394 3395 Input Parameters: 3396 . ts - The TS context obtained from TSCreate() 3397 3398 Notes: 3399 TSPostEvaluate() is typically used within time stepping implementations, 3400 most users would not generally call this routine themselves. 3401 3402 Level: developer 3403 3404 .seealso: TSSetPostEvaluate(), TSSetPreStep(), TSPreStep(), TSPostStep() 3405 @*/ 3406 PetscErrorCode TSPostEvaluate(TS ts) 3407 { 3408 PetscErrorCode ierr; 3409 3410 PetscFunctionBegin; 3411 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3412 if (ts->postevaluate) { 3413 Vec U; 3414 PetscObjectState sprev,spost; 3415 3416 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 3417 ierr = PetscObjectStateGet((PetscObject)U,&sprev);CHKERRQ(ierr); 3418 PetscStackCallStandard((*ts->postevaluate),(ts)); 3419 ierr = PetscObjectStateGet((PetscObject)U,&spost);CHKERRQ(ierr); 3420 if (sprev != spost) {ierr = TSRestartStep(ts);CHKERRQ(ierr);} 3421 } 3422 PetscFunctionReturn(0); 3423 } 3424 3425 /*@C 3426 TSSetPostStep - Sets the general-purpose function 3427 called once at the end of each time step. 3428 3429 Logically Collective on TS 3430 3431 Input Parameters: 3432 + ts - The TS context obtained from TSCreate() 3433 - func - The function 3434 3435 Calling sequence of func: 3436 $ func (TS ts); 3437 3438 Notes: 3439 The function set by TSSetPostStep() is called after each successful step. The solution vector X 3440 obtained by TSGetSolution() may be different than that computed at the step end if the event handler 3441 locates an event and TSPostEvent() modifies it. Use TSSetPostEvaluate() if an unmodified solution is needed instead. 3442 3443 Level: intermediate 3444 3445 .seealso: TSSetPreStep(), TSSetPreStage(), TSSetPostEvaluate(), TSGetTimeStep(), TSGetStepNumber(), TSGetTime(), TSRestartStep() 3446 @*/ 3447 PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS)) 3448 { 3449 PetscFunctionBegin; 3450 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 3451 ts->poststep = func; 3452 PetscFunctionReturn(0); 3453 } 3454 3455 /*@ 3456 TSPostStep - Runs the user-defined post-step function. 3457 3458 Collective on TS 3459 3460 Input Parameters: 3461 . ts - The TS context obtained from TSCreate() 3462 3463 Notes: 3464 TSPostStep() is typically used within time stepping implementations, 3465 so most users would not generally call this routine themselves. 3466 3467 Level: developer 3468 3469 @*/ 3470 PetscErrorCode TSPostStep(TS ts) 3471 { 3472 PetscErrorCode ierr; 3473 3474 PetscFunctionBegin; 3475 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3476 if (ts->poststep) { 3477 Vec U; 3478 PetscObjectState sprev,spost; 3479 3480 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 3481 ierr = PetscObjectStateGet((PetscObject)U,&sprev);CHKERRQ(ierr); 3482 PetscStackCallStandard((*ts->poststep),(ts)); 3483 ierr = PetscObjectStateGet((PetscObject)U,&spost);CHKERRQ(ierr); 3484 if (sprev != spost) {ierr = TSRestartStep(ts);CHKERRQ(ierr);} 3485 } 3486 PetscFunctionReturn(0); 3487 } 3488 3489 /*@ 3490 TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval 3491 3492 Collective on TS 3493 3494 Input Argument: 3495 + ts - time stepping context 3496 - t - time to interpolate to 3497 3498 Output Argument: 3499 . U - state at given time 3500 3501 Level: intermediate 3502 3503 Developer Notes: 3504 TSInterpolate() and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints. 3505 3506 .seealso: TSSetExactFinalTime(), TSSolve() 3507 @*/ 3508 PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec U) 3509 { 3510 PetscErrorCode ierr; 3511 3512 PetscFunctionBegin; 3513 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3514 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 3515 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); 3516 if (!ts->ops->interpolate) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name); 3517 ierr = (*ts->ops->interpolate)(ts,t,U);CHKERRQ(ierr); 3518 PetscFunctionReturn(0); 3519 } 3520 3521 /*@ 3522 TSStep - Steps one time step 3523 3524 Collective on TS 3525 3526 Input Parameter: 3527 . ts - the TS context obtained from TSCreate() 3528 3529 Level: developer 3530 3531 Notes: 3532 The public interface for the ODE/DAE solvers is TSSolve(), you should almost for sure be using that routine and not this routine. 3533 3534 The hook set using TSSetPreStep() is called before each attempt to take the step. In general, the time step size may 3535 be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages. 3536 3537 This may over-step the final time provided in TSSetMaxTime() depending on the time-step used. TSSolve() interpolates to exactly the 3538 time provided in TSSetMaxTime(). One can use TSInterpolate() to determine an interpolated solution within the final timestep. 3539 3540 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSInterpolate() 3541 @*/ 3542 PetscErrorCode TSStep(TS ts) 3543 { 3544 PetscErrorCode ierr; 3545 static PetscBool cite = PETSC_FALSE; 3546 PetscReal ptime; 3547 3548 PetscFunctionBegin; 3549 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3550 ierr = PetscCitationsRegister("@article{tspaper,\n" 3551 " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n" 3552 " author = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n" 3553 " journal = {arXiv e-preprints},\n" 3554 " eprint = {1806.01437},\n" 3555 " archivePrefix = {arXiv},\n" 3556 " year = {2018}\n}\n",&cite);CHKERRQ(ierr); 3557 3558 ierr = TSSetUp(ts);CHKERRQ(ierr); 3559 ierr = TSTrajectorySetUp(ts->trajectory,ts);CHKERRQ(ierr); 3560 3561 if (!ts->ops->step) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSStep not implemented for type '%s'",((PetscObject)ts)->type_name); 3562 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>"); 3563 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()"); 3564 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"); 3565 3566 if (!ts->steps) ts->ptime_prev = ts->ptime; 3567 ptime = ts->ptime; ts->ptime_prev_rollback = ts->ptime_prev; 3568 ts->reason = TS_CONVERGED_ITERATING; 3569 3570 ierr = PetscLogEventBegin(TS_Step,ts,0,0,0);CHKERRQ(ierr); 3571 ierr = (*ts->ops->step)(ts);CHKERRQ(ierr); 3572 ierr = PetscLogEventEnd(TS_Step,ts,0,0,0);CHKERRQ(ierr); 3573 3574 if (ts->reason >= 0) { 3575 ts->ptime_prev = ptime; 3576 ts->steps++; 3577 ts->steprollback = PETSC_FALSE; 3578 ts->steprestart = PETSC_FALSE; 3579 } 3580 3581 if (!ts->reason) { 3582 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 3583 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 3584 } 3585 3586 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]); 3587 if (ts->reason < 0 && ts->errorifstepfailed) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]); 3588 PetscFunctionReturn(0); 3589 } 3590 3591 /*@ 3592 TSEvaluateWLTE - Evaluate the weighted local truncation error norm 3593 at the end of a time step with a given order of accuracy. 3594 3595 Collective on TS 3596 3597 Input Arguments: 3598 + ts - time stepping context 3599 . wnormtype - norm type, either NORM_2 or NORM_INFINITY 3600 - order - optional, desired order for the error evaluation or PETSC_DECIDE 3601 3602 Output Arguments: 3603 + order - optional, the actual order of the error evaluation 3604 - wlte - the weighted local truncation error norm 3605 3606 Level: advanced 3607 3608 Notes: 3609 If the timestepper cannot evaluate the error in a particular step 3610 (eg. in the first step or restart steps after event handling), 3611 this routine returns wlte=-1.0 . 3612 3613 .seealso: TSStep(), TSAdapt, TSErrorWeightedNorm() 3614 @*/ 3615 PetscErrorCode TSEvaluateWLTE(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte) 3616 { 3617 PetscErrorCode ierr; 3618 3619 PetscFunctionBegin; 3620 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3621 PetscValidType(ts,1); 3622 PetscValidLogicalCollectiveEnum(ts,wnormtype,2); 3623 if (order) PetscValidIntPointer(order,3); 3624 if (order) PetscValidLogicalCollectiveInt(ts,*order,3); 3625 PetscValidRealPointer(wlte,4); 3626 if (wnormtype != NORM_2 && wnormtype != NORM_INFINITY) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"No support for norm type %s",NormTypes[wnormtype]); 3627 if (!ts->ops->evaluatewlte) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSEvaluateWLTE not implemented for type '%s'",((PetscObject)ts)->type_name); 3628 ierr = (*ts->ops->evaluatewlte)(ts,wnormtype,order,wlte);CHKERRQ(ierr); 3629 PetscFunctionReturn(0); 3630 } 3631 3632 /*@ 3633 TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy. 3634 3635 Collective on TS 3636 3637 Input Arguments: 3638 + ts - time stepping context 3639 . order - desired order of accuracy 3640 - done - whether the step was evaluated at this order (pass NULL to generate an error if not available) 3641 3642 Output Arguments: 3643 . U - state at the end of the current step 3644 3645 Level: advanced 3646 3647 Notes: 3648 This function cannot be called until all stages have been evaluated. 3649 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. 3650 3651 .seealso: TSStep(), TSAdapt 3652 @*/ 3653 PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec U,PetscBool *done) 3654 { 3655 PetscErrorCode ierr; 3656 3657 PetscFunctionBegin; 3658 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3659 PetscValidType(ts,1); 3660 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 3661 if (!ts->ops->evaluatestep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSEvaluateStep not implemented for type '%s'",((PetscObject)ts)->type_name); 3662 ierr = (*ts->ops->evaluatestep)(ts,order,U,done);CHKERRQ(ierr); 3663 PetscFunctionReturn(0); 3664 } 3665 3666 /*@C 3667 TSGetComputeInitialCondition - Get the function used to automatically compute an initial condition for the timestepping. 3668 3669 Not collective 3670 3671 Input Argument: 3672 . ts - time stepping context 3673 3674 Output Argument: 3675 . initConditions - The function which computes an initial condition 3676 3677 Level: advanced 3678 3679 Notes: 3680 The calling sequence for the function is 3681 $ initCondition(TS ts, Vec u) 3682 $ ts - The timestepping context 3683 $ u - The input vector in which the initial condition is stored 3684 3685 .seealso: TSSetComputeInitialCondition(), TSComputeInitialCondition() 3686 @*/ 3687 PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS, Vec)) 3688 { 3689 PetscFunctionBegin; 3690 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3691 PetscValidPointer(initCondition, 2); 3692 *initCondition = ts->ops->initcondition; 3693 PetscFunctionReturn(0); 3694 } 3695 3696 /*@C 3697 TSSetComputeInitialCondition - Set the function used to automatically compute an initial condition for the timestepping. 3698 3699 Logically collective on ts 3700 3701 Input Arguments: 3702 + ts - time stepping context 3703 - initCondition - The function which computes an initial condition 3704 3705 Level: advanced 3706 3707 Calling sequence for initCondition: 3708 $ PetscErrorCode initCondition(TS ts, Vec u) 3709 3710 + ts - The timestepping context 3711 - u - The input vector in which the initial condition is to be stored 3712 3713 .seealso: TSGetComputeInitialCondition(), TSComputeInitialCondition() 3714 @*/ 3715 PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS, Vec)) 3716 { 3717 PetscFunctionBegin; 3718 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3719 PetscValidFunction(initCondition, 2); 3720 ts->ops->initcondition = initCondition; 3721 PetscFunctionReturn(0); 3722 } 3723 3724 /*@ 3725 TSComputeInitialCondition - Compute an initial condition for the timestepping using the function previously set. 3726 3727 Collective on ts 3728 3729 Input Arguments: 3730 + ts - time stepping context 3731 - u - The Vec to store the condition in which will be used in TSSolve() 3732 3733 Level: advanced 3734 3735 .seealso: TSGetComputeInitialCondition(), TSSetComputeInitialCondition(), TSSolve() 3736 @*/ 3737 PetscErrorCode TSComputeInitialCondition(TS ts, Vec u) 3738 { 3739 PetscErrorCode ierr; 3740 3741 PetscFunctionBegin; 3742 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3743 PetscValidHeaderSpecific(u, VEC_CLASSID, 2); 3744 if (ts->ops->initcondition) {ierr = (*ts->ops->initcondition)(ts, u);CHKERRQ(ierr);} 3745 PetscFunctionReturn(0); 3746 } 3747 3748 /*@C 3749 TSGetComputeExactError - Get the function used to automatically compute the exact error for the timestepping. 3750 3751 Not collective 3752 3753 Input Argument: 3754 . ts - time stepping context 3755 3756 Output Argument: 3757 . exactError - The function which computes the solution error 3758 3759 Level: advanced 3760 3761 Calling sequence for exactError: 3762 $ PetscErrorCode exactError(TS ts, Vec u) 3763 3764 + ts - The timestepping context 3765 . u - The approximate solution vector 3766 - e - The input vector in which the error is stored 3767 3768 .seealso: TSGetComputeExactError(), TSComputeExactError() 3769 @*/ 3770 PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS, Vec, Vec)) 3771 { 3772 PetscFunctionBegin; 3773 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3774 PetscValidPointer(exactError, 2); 3775 *exactError = ts->ops->exacterror; 3776 PetscFunctionReturn(0); 3777 } 3778 3779 /*@C 3780 TSSetComputeExactError - Set the function used to automatically compute the exact error for the timestepping. 3781 3782 Logically collective on ts 3783 3784 Input Arguments: 3785 + ts - time stepping context 3786 - exactError - The function which computes the solution error 3787 3788 Level: advanced 3789 3790 Calling sequence for exactError: 3791 $ PetscErrorCode exactError(TS ts, Vec u) 3792 3793 + ts - The timestepping context 3794 . u - The approximate solution vector 3795 - e - The input vector in which the error is stored 3796 3797 .seealso: TSGetComputeExactError(), TSComputeExactError() 3798 @*/ 3799 PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS, Vec, Vec)) 3800 { 3801 PetscFunctionBegin; 3802 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3803 PetscValidFunction(exactError, 2); 3804 ts->ops->exacterror = exactError; 3805 PetscFunctionReturn(0); 3806 } 3807 3808 /*@ 3809 TSComputeExactError - Compute the solution error for the timestepping using the function previously set. 3810 3811 Collective on ts 3812 3813 Input Arguments: 3814 + ts - time stepping context 3815 . u - The approximate solution 3816 - e - The Vec used to store the error 3817 3818 Level: advanced 3819 3820 .seealso: TSGetComputeInitialCondition(), TSSetComputeInitialCondition(), TSSolve() 3821 @*/ 3822 PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e) 3823 { 3824 PetscErrorCode ierr; 3825 3826 PetscFunctionBegin; 3827 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3828 PetscValidHeaderSpecific(u, VEC_CLASSID, 2); 3829 PetscValidHeaderSpecific(e, VEC_CLASSID, 3); 3830 if (ts->ops->exacterror) {ierr = (*ts->ops->exacterror)(ts, u, e);CHKERRQ(ierr);} 3831 PetscFunctionReturn(0); 3832 } 3833 3834 /*@ 3835 TSSolve - Steps the requested number of timesteps. 3836 3837 Collective on TS 3838 3839 Input Parameter: 3840 + ts - the TS context obtained from TSCreate() 3841 - u - the solution vector (can be null if TSSetSolution() was used and TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP) was not used, 3842 otherwise must contain the initial conditions and will contain the solution at the final requested time 3843 3844 Level: beginner 3845 3846 Notes: 3847 The final time returned by this function may be different from the time of the internally 3848 held state accessible by TSGetSolution() and TSGetTime() because the method may have 3849 stepped over the final time. 3850 3851 .seealso: TSCreate(), TSSetSolution(), TSStep(), TSGetTime(), TSGetSolveTime() 3852 @*/ 3853 PetscErrorCode TSSolve(TS ts,Vec u) 3854 { 3855 Vec solution; 3856 PetscErrorCode ierr; 3857 3858 PetscFunctionBegin; 3859 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3860 if (u) PetscValidHeaderSpecific(u,VEC_CLASSID,2); 3861 3862 ierr = TSSetExactFinalTimeDefault(ts);CHKERRQ(ierr); 3863 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 */ 3864 if (!ts->vec_sol || u == ts->vec_sol) { 3865 ierr = VecDuplicate(u,&solution);CHKERRQ(ierr); 3866 ierr = TSSetSolution(ts,solution);CHKERRQ(ierr); 3867 ierr = VecDestroy(&solution);CHKERRQ(ierr); /* grant ownership */ 3868 } 3869 ierr = VecCopy(u,ts->vec_sol);CHKERRQ(ierr); 3870 if (ts->forward_solve) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE"); 3871 } else if (u) { 3872 ierr = TSSetSolution(ts,u);CHKERRQ(ierr); 3873 } 3874 ierr = TSSetUp(ts);CHKERRQ(ierr); 3875 ierr = TSTrajectorySetUp(ts->trajectory,ts);CHKERRQ(ierr); 3876 3877 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>"); 3878 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()"); 3879 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"); 3880 3881 if (ts->forward_solve) { 3882 ierr = TSForwardSetUp(ts);CHKERRQ(ierr); 3883 } 3884 3885 /* reset number of steps only when the step is not restarted. ARKIMEX 3886 restarts the step after an event. Resetting these counters in such case causes 3887 TSTrajectory to incorrectly save the output files 3888 */ 3889 /* reset time step and iteration counters */ 3890 if (!ts->steps) { 3891 ts->ksp_its = 0; 3892 ts->snes_its = 0; 3893 ts->num_snes_failures = 0; 3894 ts->reject = 0; 3895 ts->steprestart = PETSC_TRUE; 3896 ts->steprollback = PETSC_FALSE; 3897 ts->rhsjacobian.time = PETSC_MIN_REAL; 3898 } 3899 3900 /* make sure initial time step does not overshoot final time */ 3901 if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) { 3902 PetscReal maxdt = ts->max_time-ts->ptime; 3903 PetscReal dt = ts->time_step; 3904 3905 ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt,maxdt,10*PETSC_MACHINE_EPSILON,0) ? maxdt : dt); 3906 } 3907 ts->reason = TS_CONVERGED_ITERATING; 3908 3909 { 3910 PetscViewer viewer; 3911 PetscViewerFormat format; 3912 PetscBool flg; 3913 static PetscBool incall = PETSC_FALSE; 3914 3915 if (!incall) { 3916 /* Estimate the convergence rate of the time discretization */ 3917 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject) ts),((PetscObject)ts)->options, ((PetscObject) ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg);CHKERRQ(ierr); 3918 if (flg) { 3919 PetscConvEst conv; 3920 DM dm; 3921 PetscReal *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */ 3922 PetscInt Nf; 3923 PetscBool checkTemporal = PETSC_TRUE; 3924 3925 incall = PETSC_TRUE; 3926 ierr = PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject) ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg);CHKERRQ(ierr); 3927 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 3928 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 3929 ierr = PetscCalloc1(PetscMax(Nf, 1), &alpha);CHKERRQ(ierr); 3930 ierr = PetscConvEstCreate(PetscObjectComm((PetscObject) ts), &conv);CHKERRQ(ierr); 3931 ierr = PetscConvEstUseTS(conv, checkTemporal);CHKERRQ(ierr); 3932 ierr = PetscConvEstSetSolver(conv, (PetscObject) ts);CHKERRQ(ierr); 3933 ierr = PetscConvEstSetFromOptions(conv);CHKERRQ(ierr); 3934 ierr = PetscConvEstSetUp(conv);CHKERRQ(ierr); 3935 ierr = PetscConvEstGetConvRate(conv, alpha);CHKERRQ(ierr); 3936 ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr); 3937 ierr = PetscConvEstRateView(conv, alpha, viewer);CHKERRQ(ierr); 3938 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 3939 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3940 ierr = PetscConvEstDestroy(&conv);CHKERRQ(ierr); 3941 ierr = PetscFree(alpha);CHKERRQ(ierr); 3942 incall = PETSC_FALSE; 3943 } 3944 } 3945 } 3946 3947 ierr = TSViewFromOptions(ts,NULL,"-ts_view_pre");CHKERRQ(ierr); 3948 3949 if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */ 3950 ierr = (*ts->ops->solve)(ts);CHKERRQ(ierr); 3951 if (u) {ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr);} 3952 ts->solvetime = ts->ptime; 3953 solution = ts->vec_sol; 3954 } else { /* Step the requested number of timesteps. */ 3955 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 3956 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 3957 3958 if (!ts->steps) { 3959 ierr = TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 3960 ierr = TSEventInitialize(ts->event,ts,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 3961 } 3962 3963 while (!ts->reason) { 3964 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 3965 if (!ts->steprollback) { 3966 ierr = TSPreStep(ts);CHKERRQ(ierr); 3967 } 3968 ierr = TSStep(ts);CHKERRQ(ierr); 3969 if (ts->testjacobian) { 3970 ierr = TSRHSJacobianTest(ts,NULL);CHKERRQ(ierr); 3971 } 3972 if (ts->testjacobiantranspose) { 3973 ierr = TSRHSJacobianTestTranspose(ts,NULL);CHKERRQ(ierr); 3974 } 3975 if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */ 3976 if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */ 3977 ierr = TSForwardCostIntegral(ts);CHKERRQ(ierr); 3978 if (ts->reason >= 0) ts->steps++; 3979 } 3980 if (ts->forward_solve) { /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */ 3981 if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */ 3982 ierr = TSForwardStep(ts);CHKERRQ(ierr); 3983 if (ts->reason >= 0) ts->steps++; 3984 } 3985 ierr = TSPostEvaluate(ts);CHKERRQ(ierr); 3986 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. */ 3987 if (ts->steprollback) { 3988 ierr = TSPostEvaluate(ts);CHKERRQ(ierr); 3989 } 3990 if (!ts->steprollback) { 3991 ierr = TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 3992 ierr = TSPostStep(ts);CHKERRQ(ierr); 3993 } 3994 } 3995 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 3996 3997 if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) { 3998 ierr = TSInterpolate(ts,ts->max_time,u);CHKERRQ(ierr); 3999 ts->solvetime = ts->max_time; 4000 solution = u; 4001 ierr = TSMonitor(ts,-1,ts->solvetime,solution);CHKERRQ(ierr); 4002 } else { 4003 if (u) {ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr);} 4004 ts->solvetime = ts->ptime; 4005 solution = ts->vec_sol; 4006 } 4007 } 4008 4009 ierr = TSViewFromOptions(ts,NULL,"-ts_view");CHKERRQ(ierr); 4010 ierr = VecViewFromOptions(solution,(PetscObject)ts,"-ts_view_solution");CHKERRQ(ierr); 4011 ierr = PetscObjectSAWsBlock((PetscObject)ts);CHKERRQ(ierr); 4012 if (ts->adjoint_solve) { 4013 ierr = TSAdjointSolve(ts);CHKERRQ(ierr); 4014 } 4015 PetscFunctionReturn(0); 4016 } 4017 4018 /*@ 4019 TSGetTime - Gets the time of the most recently completed step. 4020 4021 Not Collective 4022 4023 Input Parameter: 4024 . ts - the TS context obtained from TSCreate() 4025 4026 Output Parameter: 4027 . t - the current time. This time may not corresponds to the final time set with TSSetMaxTime(), use TSGetSolveTime(). 4028 4029 Level: beginner 4030 4031 Note: 4032 When called during time step evaluation (e.g. during residual evaluation or via hooks set using TSSetPreStep(), 4033 TSSetPreStage(), TSSetPostStage(), or TSSetPostStep()), the time is the time at the start of the step being evaluated. 4034 4035 .seealso: TSGetSolveTime(), TSSetTime(), TSGetTimeStep(), TSGetStepNumber() 4036 4037 @*/ 4038 PetscErrorCode TSGetTime(TS ts,PetscReal *t) 4039 { 4040 PetscFunctionBegin; 4041 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4042 PetscValidRealPointer(t,2); 4043 *t = ts->ptime; 4044 PetscFunctionReturn(0); 4045 } 4046 4047 /*@ 4048 TSGetPrevTime - Gets the starting time of the previously completed step. 4049 4050 Not Collective 4051 4052 Input Parameter: 4053 . ts - the TS context obtained from TSCreate() 4054 4055 Output Parameter: 4056 . t - the previous time 4057 4058 Level: beginner 4059 4060 .seealso: TSGetTime(), TSGetSolveTime(), TSGetTimeStep() 4061 4062 @*/ 4063 PetscErrorCode TSGetPrevTime(TS ts,PetscReal *t) 4064 { 4065 PetscFunctionBegin; 4066 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4067 PetscValidRealPointer(t,2); 4068 *t = ts->ptime_prev; 4069 PetscFunctionReturn(0); 4070 } 4071 4072 /*@ 4073 TSSetTime - Allows one to reset the time. 4074 4075 Logically Collective on TS 4076 4077 Input Parameters: 4078 + ts - the TS context obtained from TSCreate() 4079 - time - the time 4080 4081 Level: intermediate 4082 4083 .seealso: TSGetTime(), TSSetMaxSteps() 4084 4085 @*/ 4086 PetscErrorCode TSSetTime(TS ts, PetscReal t) 4087 { 4088 PetscFunctionBegin; 4089 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4090 PetscValidLogicalCollectiveReal(ts,t,2); 4091 ts->ptime = t; 4092 PetscFunctionReturn(0); 4093 } 4094 4095 /*@C 4096 TSSetOptionsPrefix - Sets the prefix used for searching for all 4097 TS options in the database. 4098 4099 Logically Collective on TS 4100 4101 Input Parameter: 4102 + ts - The TS context 4103 - prefix - The prefix to prepend to all option names 4104 4105 Notes: 4106 A hyphen (-) must NOT be given at the beginning of the prefix name. 4107 The first character of all runtime options is AUTOMATICALLY the 4108 hyphen. 4109 4110 Level: advanced 4111 4112 .seealso: TSSetFromOptions() 4113 4114 @*/ 4115 PetscErrorCode TSSetOptionsPrefix(TS ts,const char prefix[]) 4116 { 4117 PetscErrorCode ierr; 4118 SNES snes; 4119 4120 PetscFunctionBegin; 4121 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4122 ierr = PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 4123 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 4124 ierr = SNESSetOptionsPrefix(snes,prefix);CHKERRQ(ierr); 4125 PetscFunctionReturn(0); 4126 } 4127 4128 /*@C 4129 TSAppendOptionsPrefix - Appends to the prefix used for searching for all 4130 TS options in the database. 4131 4132 Logically Collective on TS 4133 4134 Input Parameter: 4135 + ts - The TS context 4136 - prefix - The prefix to prepend to all option names 4137 4138 Notes: 4139 A hyphen (-) must NOT be given at the beginning of the prefix name. 4140 The first character of all runtime options is AUTOMATICALLY the 4141 hyphen. 4142 4143 Level: advanced 4144 4145 .seealso: TSGetOptionsPrefix() 4146 4147 @*/ 4148 PetscErrorCode TSAppendOptionsPrefix(TS ts,const char prefix[]) 4149 { 4150 PetscErrorCode ierr; 4151 SNES snes; 4152 4153 PetscFunctionBegin; 4154 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4155 ierr = PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 4156 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 4157 ierr = SNESAppendOptionsPrefix(snes,prefix);CHKERRQ(ierr); 4158 PetscFunctionReturn(0); 4159 } 4160 4161 /*@C 4162 TSGetOptionsPrefix - Sets the prefix used for searching for all 4163 TS options in the database. 4164 4165 Not Collective 4166 4167 Input Parameter: 4168 . ts - The TS context 4169 4170 Output Parameter: 4171 . prefix - A pointer to the prefix string used 4172 4173 Notes: 4174 On the fortran side, the user should pass in a string 'prifix' of 4175 sufficient length to hold the prefix. 4176 4177 Level: intermediate 4178 4179 .seealso: TSAppendOptionsPrefix() 4180 @*/ 4181 PetscErrorCode TSGetOptionsPrefix(TS ts,const char *prefix[]) 4182 { 4183 PetscErrorCode ierr; 4184 4185 PetscFunctionBegin; 4186 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4187 PetscValidPointer(prefix,2); 4188 ierr = PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 4189 PetscFunctionReturn(0); 4190 } 4191 4192 /*@C 4193 TSGetRHSJacobian - Returns the Jacobian J at the present timestep. 4194 4195 Not Collective, but parallel objects are returned if TS is parallel 4196 4197 Input Parameter: 4198 . ts - The TS context obtained from TSCreate() 4199 4200 Output Parameters: 4201 + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or NULL) 4202 . Pmat - The matrix from which the preconditioner is constructed, usually the same as Amat (or NULL) 4203 . func - Function to compute the Jacobian of the RHS (or NULL) 4204 - ctx - User-defined context for Jacobian evaluation routine (or NULL) 4205 4206 Notes: 4207 You can pass in NULL for any return argument you do not need. 4208 4209 Level: intermediate 4210 4211 .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetStepNumber() 4212 4213 @*/ 4214 PetscErrorCode TSGetRHSJacobian(TS ts,Mat *Amat,Mat *Pmat,TSRHSJacobian *func,void **ctx) 4215 { 4216 PetscErrorCode ierr; 4217 DM dm; 4218 4219 PetscFunctionBegin; 4220 if (Amat || Pmat) { 4221 SNES snes; 4222 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 4223 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 4224 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 4225 } 4226 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 4227 ierr = DMTSGetRHSJacobian(dm,func,ctx);CHKERRQ(ierr); 4228 PetscFunctionReturn(0); 4229 } 4230 4231 /*@C 4232 TSGetIJacobian - Returns the implicit Jacobian at the present timestep. 4233 4234 Not Collective, but parallel objects are returned if TS is parallel 4235 4236 Input Parameter: 4237 . ts - The TS context obtained from TSCreate() 4238 4239 Output Parameters: 4240 + Amat - The (approximate) Jacobian of F(t,U,U_t) 4241 . Pmat - The matrix from which the preconditioner is constructed, often the same as Amat 4242 . f - The function to compute the matrices 4243 - ctx - User-defined context for Jacobian evaluation routine 4244 4245 Notes: 4246 You can pass in NULL for any return argument you do not need. 4247 4248 Level: advanced 4249 4250 .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetStepNumber() 4251 4252 @*/ 4253 PetscErrorCode TSGetIJacobian(TS ts,Mat *Amat,Mat *Pmat,TSIJacobian *f,void **ctx) 4254 { 4255 PetscErrorCode ierr; 4256 DM dm; 4257 4258 PetscFunctionBegin; 4259 if (Amat || Pmat) { 4260 SNES snes; 4261 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 4262 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 4263 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 4264 } 4265 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 4266 ierr = DMTSGetIJacobian(dm,f,ctx);CHKERRQ(ierr); 4267 PetscFunctionReturn(0); 4268 } 4269 4270 #include <petsc/private/dmimpl.h> 4271 /*@ 4272 TSSetDM - Sets the DM that may be used by some nonlinear solvers or preconditioners under the TS 4273 4274 Logically Collective on ts 4275 4276 Input Parameters: 4277 + ts - the ODE integrator object 4278 - dm - the dm, cannot be NULL 4279 4280 Notes: 4281 A DM can only be used for solving one problem at a time because information about the problem is stored on the DM, 4282 even when not using interfaces like DMTSSetIFunction(). Use DMClone() to get a distinct DM when solving 4283 different problems using the same function space. 4284 4285 Level: intermediate 4286 4287 .seealso: TSGetDM(), SNESSetDM(), SNESGetDM() 4288 @*/ 4289 PetscErrorCode TSSetDM(TS ts,DM dm) 4290 { 4291 PetscErrorCode ierr; 4292 SNES snes; 4293 DMTS tsdm; 4294 4295 PetscFunctionBegin; 4296 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4297 PetscValidHeaderSpecific(dm,DM_CLASSID,2); 4298 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 4299 if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */ 4300 if (ts->dm->dmts && !dm->dmts) { 4301 ierr = DMCopyDMTS(ts->dm,dm);CHKERRQ(ierr); 4302 ierr = DMGetDMTS(ts->dm,&tsdm);CHKERRQ(ierr); 4303 if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */ 4304 tsdm->originaldm = dm; 4305 } 4306 } 4307 ierr = DMDestroy(&ts->dm);CHKERRQ(ierr); 4308 } 4309 ts->dm = dm; 4310 4311 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 4312 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 4313 PetscFunctionReturn(0); 4314 } 4315 4316 /*@ 4317 TSGetDM - Gets the DM that may be used by some preconditioners 4318 4319 Not Collective 4320 4321 Input Parameter: 4322 . ts - the preconditioner context 4323 4324 Output Parameter: 4325 . dm - the dm 4326 4327 Level: intermediate 4328 4329 .seealso: TSSetDM(), SNESSetDM(), SNESGetDM() 4330 @*/ 4331 PetscErrorCode TSGetDM(TS ts,DM *dm) 4332 { 4333 PetscErrorCode ierr; 4334 4335 PetscFunctionBegin; 4336 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4337 if (!ts->dm) { 4338 ierr = DMShellCreate(PetscObjectComm((PetscObject)ts),&ts->dm);CHKERRQ(ierr); 4339 if (ts->snes) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 4340 } 4341 *dm = ts->dm; 4342 PetscFunctionReturn(0); 4343 } 4344 4345 /*@ 4346 SNESTSFormFunction - Function to evaluate nonlinear residual 4347 4348 Logically Collective on SNES 4349 4350 Input Parameter: 4351 + snes - nonlinear solver 4352 . U - the current state at which to evaluate the residual 4353 - ctx - user context, must be a TS 4354 4355 Output Parameter: 4356 . F - the nonlinear residual 4357 4358 Notes: 4359 This function is not normally called by users and is automatically registered with the SNES used by TS. 4360 It is most frequently passed to MatFDColoringSetFunction(). 4361 4362 Level: advanced 4363 4364 .seealso: SNESSetFunction(), MatFDColoringSetFunction() 4365 @*/ 4366 PetscErrorCode SNESTSFormFunction(SNES snes,Vec U,Vec F,void *ctx) 4367 { 4368 TS ts = (TS)ctx; 4369 PetscErrorCode ierr; 4370 4371 PetscFunctionBegin; 4372 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4373 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 4374 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 4375 PetscValidHeaderSpecific(ts,TS_CLASSID,4); 4376 ierr = (ts->ops->snesfunction)(snes,U,F,ts);CHKERRQ(ierr); 4377 PetscFunctionReturn(0); 4378 } 4379 4380 /*@ 4381 SNESTSFormJacobian - Function to evaluate the Jacobian 4382 4383 Collective on SNES 4384 4385 Input Parameter: 4386 + snes - nonlinear solver 4387 . U - the current state at which to evaluate the residual 4388 - ctx - user context, must be a TS 4389 4390 Output Parameter: 4391 + A - the Jacobian 4392 . B - the preconditioning matrix (may be the same as A) 4393 - flag - indicates any structure change in the matrix 4394 4395 Notes: 4396 This function is not normally called by users and is automatically registered with the SNES used by TS. 4397 4398 Level: developer 4399 4400 .seealso: SNESSetJacobian() 4401 @*/ 4402 PetscErrorCode SNESTSFormJacobian(SNES snes,Vec U,Mat A,Mat B,void *ctx) 4403 { 4404 TS ts = (TS)ctx; 4405 PetscErrorCode ierr; 4406 4407 PetscFunctionBegin; 4408 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4409 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 4410 PetscValidPointer(A,3); 4411 PetscValidHeaderSpecific(A,MAT_CLASSID,3); 4412 PetscValidPointer(B,4); 4413 PetscValidHeaderSpecific(B,MAT_CLASSID,4); 4414 PetscValidHeaderSpecific(ts,TS_CLASSID,5); 4415 ierr = (ts->ops->snesjacobian)(snes,U,A,B,ts);CHKERRQ(ierr); 4416 PetscFunctionReturn(0); 4417 } 4418 4419 /*@C 4420 TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems Udot = A U only 4421 4422 Collective on TS 4423 4424 Input Arguments: 4425 + ts - time stepping context 4426 . t - time at which to evaluate 4427 . U - state at which to evaluate 4428 - ctx - context 4429 4430 Output Arguments: 4431 . F - right hand side 4432 4433 Level: intermediate 4434 4435 Notes: 4436 This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems. 4437 The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian(). 4438 4439 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant() 4440 @*/ 4441 PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec U,Vec F,void *ctx) 4442 { 4443 PetscErrorCode ierr; 4444 Mat Arhs,Brhs; 4445 4446 PetscFunctionBegin; 4447 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 4448 /* undo the damage caused by shifting */ 4449 ierr = TSRecoverRHSJacobian(ts,Arhs,Brhs);CHKERRQ(ierr); 4450 ierr = TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);CHKERRQ(ierr); 4451 ierr = MatMult(Arhs,U,F);CHKERRQ(ierr); 4452 PetscFunctionReturn(0); 4453 } 4454 4455 /*@C 4456 TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent. 4457 4458 Collective on TS 4459 4460 Input Arguments: 4461 + ts - time stepping context 4462 . t - time at which to evaluate 4463 . U - state at which to evaluate 4464 - ctx - context 4465 4466 Output Arguments: 4467 + A - pointer to operator 4468 . B - pointer to preconditioning matrix 4469 - flg - matrix structure flag 4470 4471 Level: intermediate 4472 4473 Notes: 4474 This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems. 4475 4476 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear() 4477 @*/ 4478 PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx) 4479 { 4480 PetscFunctionBegin; 4481 PetscFunctionReturn(0); 4482 } 4483 4484 /*@C 4485 TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only 4486 4487 Collective on TS 4488 4489 Input Arguments: 4490 + ts - time stepping context 4491 . t - time at which to evaluate 4492 . U - state at which to evaluate 4493 . Udot - time derivative of state vector 4494 - ctx - context 4495 4496 Output Arguments: 4497 . F - left hand side 4498 4499 Level: intermediate 4500 4501 Notes: 4502 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 4503 user is required to write their own TSComputeIFunction. 4504 This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems. 4505 The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian(). 4506 4507 Note that using this function is NOT equivalent to using TSComputeRHSFunctionLinear() since that solves Udot = A U 4508 4509 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant(), TSComputeRHSFunctionLinear() 4510 @*/ 4511 PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx) 4512 { 4513 PetscErrorCode ierr; 4514 Mat A,B; 4515 4516 PetscFunctionBegin; 4517 ierr = TSGetIJacobian(ts,&A,&B,NULL,NULL);CHKERRQ(ierr); 4518 ierr = TSComputeIJacobian(ts,t,U,Udot,1.0,A,B,PETSC_TRUE);CHKERRQ(ierr); 4519 ierr = MatMult(A,Udot,F);CHKERRQ(ierr); 4520 PetscFunctionReturn(0); 4521 } 4522 4523 /*@C 4524 TSComputeIJacobianConstant - Reuses a time-independent for a semi-implicit DAE or ODE 4525 4526 Collective on TS 4527 4528 Input Arguments: 4529 + ts - time stepping context 4530 . t - time at which to evaluate 4531 . U - state at which to evaluate 4532 . Udot - time derivative of state vector 4533 . shift - shift to apply 4534 - ctx - context 4535 4536 Output Arguments: 4537 + A - pointer to operator 4538 . B - pointer to preconditioning matrix 4539 - flg - matrix structure flag 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 Parameter: 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 Parameter: 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 Parameter: 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 Parameter: 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 Arguments: 4890 . ts - time stepping context 4891 4892 Output Arguments: 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 Arguments: 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 Arguments: 4969 . ts - time integration context 4970 4971 Output Arguments: 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 Arguments: 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 Arguments: 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 Arguments: 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 Arguments: 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 Arguments: 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 Arguments: 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 Arguments: 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 Arguments: 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 Arguments: 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 Arguments: 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 Arguments: 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 Arguments: 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 Arguments: 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 Arguments: 5756 . ts - time stepping context 5757 5758 Output Arguments: 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 Parameter: 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