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