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