1 /* 2 Provides a PETSc interface to SUNDIALS/CVODE solver. 3 The interface to PVODE (old version of CVODE) was originally contributed 4 by Liyang Xu. It has been redone by Hong Zhang and Dinesh Kaushik. 5 6 Reference: sundials-2.4.0/examples/cvode/parallel/cvDiurnal_kry_p.c 7 */ 8 #include "sundials.h" /*I "petscts.h" I*/ 9 10 /* 11 TSPrecond_Sundials - function that we provide to SUNDIALS to 12 evaluate the preconditioner. 13 */ 14 #undef __FUNCT__ 15 #define __FUNCT__ "TSPrecond_Sundials" 16 PetscErrorCode TSPrecond_Sundials(realtype tn,N_Vector y,N_Vector fy,booleantype jok,booleantype *jcurPtr, 17 realtype _gamma,void *P_data,N_Vector vtemp1,N_Vector vtemp2,N_Vector vtemp3) 18 { 19 TS ts = (TS) P_data; 20 TS_Sundials *cvode = (TS_Sundials*)ts->data; 21 PC pc; 22 PetscErrorCode ierr; 23 Mat J,P; 24 Vec yy = cvode->w1,yydot = cvode->ydot; 25 PetscReal gm = (PetscReal)_gamma; 26 MatStructure str = DIFFERENT_NONZERO_PATTERN; 27 PetscScalar *y_data; 28 29 PetscFunctionBegin; 30 ierr = TSGetIJacobian(ts,&J,&P,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 31 y_data = (PetscScalar*) N_VGetArrayPointer(y); 32 ierr = VecPlaceArray(yy,y_data);CHKERRQ(ierr); 33 ierr = VecZeroEntries(yydot);CHKERRQ(ierr); /* The Jacobian is independent of Ydot for ODE which is all that CVode works for */ 34 /* compute the shifted Jacobian (1/gm)*I + Jrest */ 35 ierr = TSComputeIJacobian(ts,ts->ptime,yy,yydot,1/gm,&J,&P,&str,PETSC_FALSE);CHKERRQ(ierr); 36 ierr = VecResetArray(yy);CHKERRQ(ierr); 37 ierr = MatScale(P,gm);CHKERRQ(ierr); /* turn into I-gm*Jrest, J is not used by Sundials */ 38 *jcurPtr = TRUE; 39 ierr = TSSundialsGetPC(ts,&pc);CHKERRQ(ierr); 40 ierr = PCSetOperators(pc,J,P,str);CHKERRQ(ierr); 41 PetscFunctionReturn(0); 42 } 43 44 /* 45 TSPSolve_Sundials - routine that we provide to Sundials that applies the preconditioner. 46 */ 47 #undef __FUNCT__ 48 #define __FUNCT__ "TSPSolve_Sundials" 49 PetscErrorCode TSPSolve_Sundials(realtype tn,N_Vector y,N_Vector fy,N_Vector r,N_Vector z, 50 realtype _gamma,realtype delta,int lr,void *P_data,N_Vector vtemp) 51 { 52 TS ts = (TS) P_data; 53 TS_Sundials *cvode = (TS_Sundials*)ts->data; 54 PC pc; 55 Vec rr = cvode->w1,zz = cvode->w2; 56 PetscErrorCode ierr; 57 PetscScalar *r_data,*z_data; 58 59 PetscFunctionBegin; 60 /* Make the PETSc work vectors rr and zz point to the arrays in the SUNDIALS vectors r and z respectively*/ 61 r_data = (PetscScalar*) N_VGetArrayPointer(r); 62 z_data = (PetscScalar*) N_VGetArrayPointer(z); 63 ierr = VecPlaceArray(rr,r_data);CHKERRQ(ierr); 64 ierr = VecPlaceArray(zz,z_data);CHKERRQ(ierr); 65 66 /* Solve the Px=r and put the result in zz */ 67 ierr = TSSundialsGetPC(ts,&pc);CHKERRQ(ierr); 68 ierr = PCApply(pc,rr,zz);CHKERRQ(ierr); 69 ierr = VecResetArray(rr);CHKERRQ(ierr); 70 ierr = VecResetArray(zz);CHKERRQ(ierr); 71 PetscFunctionReturn(0); 72 } 73 74 /* 75 TSFunction_Sundials - routine that we provide to Sundials that applies the right hand side. 76 */ 77 #undef __FUNCT__ 78 #define __FUNCT__ "TSFunction_Sundials" 79 int TSFunction_Sundials(realtype t,N_Vector y,N_Vector ydot,void *ctx) 80 { 81 TS ts = (TS) ctx; 82 DM dm; 83 DMTS tsdm; 84 TSIFunction ifunction; 85 MPI_Comm comm = ((PetscObject)ts)->comm; 86 TS_Sundials *cvode = (TS_Sundials*)ts->data; 87 Vec yy = cvode->w1,yyd = cvode->w2,yydot = cvode->ydot; 88 PetscScalar *y_data,*ydot_data; 89 PetscErrorCode ierr; 90 91 PetscFunctionBegin; 92 /* Make the PETSc work vectors yy and yyd point to the arrays in the SUNDIALS vectors y and ydot respectively*/ 93 y_data = (PetscScalar*) N_VGetArrayPointer(y); 94 ydot_data = (PetscScalar*) N_VGetArrayPointer(ydot); 95 ierr = VecPlaceArray(yy,y_data);CHKERRABORT(comm,ierr); 96 ierr = VecPlaceArray(yyd,ydot_data);CHKERRABORT(comm,ierr); 97 98 /* Now compute the right hand side function, via IFunction unless only the more efficient RHSFunction is set */ 99 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 100 ierr = DMGetDMTS(dm,&tsdm);CHKERRQ(ierr); 101 ierr = DMTSGetIFunction(dm,&ifunction,PETSC_NULL);CHKERRQ(ierr); 102 if (!ifunction) { 103 ierr = TSComputeRHSFunction(ts,t,yy,yyd);CHKERRQ(ierr); 104 } else { /* If rhsfunction is also set, this computes both parts and shifts them to the right */ 105 ierr = VecZeroEntries(yydot);CHKERRQ(ierr); 106 ierr = TSComputeIFunction(ts,t,yy,yydot,yyd,PETSC_FALSE);CHKERRABORT(comm,ierr); 107 ierr = VecScale(yyd,-1.);CHKERRQ(ierr); 108 } 109 ierr = VecResetArray(yy);CHKERRABORT(comm,ierr); 110 ierr = VecResetArray(yyd);CHKERRABORT(comm,ierr); 111 PetscFunctionReturn(0); 112 } 113 114 /* 115 TSStep_Sundials - Calls Sundials to integrate the ODE. 116 */ 117 #undef __FUNCT__ 118 #define __FUNCT__ "TSStep_Sundials" 119 PetscErrorCode TSStep_Sundials(TS ts) 120 { 121 TS_Sundials *cvode = (TS_Sundials*)ts->data; 122 PetscErrorCode ierr; 123 PetscInt flag; 124 long int its,nsteps; 125 realtype t,tout; 126 PetscScalar *y_data; 127 void *mem; 128 129 PetscFunctionBegin; 130 mem = cvode->mem; 131 tout = ts->max_time; 132 ierr = VecGetArray(ts->vec_sol,&y_data);CHKERRQ(ierr); 133 N_VSetArrayPointer((realtype*)y_data,cvode->y); 134 ierr = VecRestoreArray(ts->vec_sol,PETSC_NULL);CHKERRQ(ierr); 135 136 ierr = TSPreStep(ts);CHKERRQ(ierr); 137 138 /* We would like to call TSPreStep() when starting each step (including rejections) and TSPreStage() before each 139 * stage solve, but CVode does not appear to support this. */ 140 if (cvode->monitorstep) flag = CVode(mem,tout,cvode->y,&t,CV_ONE_STEP); 141 else flag = CVode(mem,tout,cvode->y,&t,CV_NORMAL); 142 143 if (flag) { /* display error message */ 144 switch (flag) { 145 case CV_ILL_INPUT: 146 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ILL_INPUT"); 147 break; 148 case CV_TOO_CLOSE: 149 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_CLOSE"); 150 break; 151 case CV_TOO_MUCH_WORK: { 152 PetscReal tcur; 153 ierr = CVodeGetNumSteps(mem,&nsteps);CHKERRQ(ierr); 154 ierr = CVodeGetCurrentTime(mem,&tcur);CHKERRQ(ierr); 155 SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_MUCH_WORK. At t=%G, nsteps %D exceeds mxstep %D. Increase '-ts_max_steps <>' or modify TSSetDuration()",tcur,nsteps,ts->max_steps); 156 } break; 157 case CV_TOO_MUCH_ACC: 158 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_MUCH_ACC"); 159 break; 160 case CV_ERR_FAILURE: 161 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ERR_FAILURE"); 162 break; 163 case CV_CONV_FAILURE: 164 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_CONV_FAILURE"); 165 break; 166 case CV_LINIT_FAIL: 167 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LINIT_FAIL"); 168 break; 169 case CV_LSETUP_FAIL: 170 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSETUP_FAIL"); 171 break; 172 case CV_LSOLVE_FAIL: 173 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSOLVE_FAIL"); 174 break; 175 case CV_RHSFUNC_FAIL: 176 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RHSFUNC_FAIL"); 177 break; 178 case CV_FIRST_RHSFUNC_ERR: 179 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_FIRST_RHSFUNC_ERR"); 180 break; 181 case CV_REPTD_RHSFUNC_ERR: 182 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_REPTD_RHSFUNC_ERR"); 183 break; 184 case CV_UNREC_RHSFUNC_ERR: 185 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_UNREC_RHSFUNC_ERR"); 186 break; 187 case CV_RTFUNC_FAIL: 188 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RTFUNC_FAIL"); 189 break; 190 default: 191 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, flag %d",flag); 192 } 193 } 194 195 /* copy the solution from cvode->y to cvode->update and sol */ 196 ierr = VecPlaceArray(cvode->w1,y_data);CHKERRQ(ierr); 197 ierr = VecCopy(cvode->w1,cvode->update);CHKERRQ(ierr); 198 ierr = VecResetArray(cvode->w1);CHKERRQ(ierr); 199 ierr = VecCopy(cvode->update,ts->vec_sol);CHKERRQ(ierr); 200 ierr = CVodeGetNumNonlinSolvIters(mem,&its);CHKERRQ(ierr); 201 ierr = CVSpilsGetNumLinIters(mem, &its); 202 ts->snes_its = its; ts->ksp_its = its; 203 204 ts->time_step = t - ts->ptime; 205 ts->ptime = t; 206 ts->steps++; 207 208 ierr = CVodeGetNumSteps(mem,&nsteps);CHKERRQ(ierr); 209 if (!cvode->monitorstep) ts->steps = nsteps; 210 PetscFunctionReturn(0); 211 } 212 213 #undef __FUNCT__ 214 #define __FUNCT__ "TSInterpolate_Sundials" 215 static PetscErrorCode TSInterpolate_Sundials(TS ts,PetscReal t,Vec X) 216 { 217 TS_Sundials *cvode = (TS_Sundials*)ts->data; 218 N_Vector y; 219 PetscErrorCode ierr; 220 PetscScalar *x_data; 221 PetscInt glosize,locsize; 222 223 PetscFunctionBegin; 224 /* get the vector size */ 225 ierr = VecGetSize(X,&glosize);CHKERRQ(ierr); 226 ierr = VecGetLocalSize(X,&locsize);CHKERRQ(ierr); 227 228 /* allocate the memory for N_Vec y */ 229 y = N_VNew_Parallel(cvode->comm_sundials,locsize,glosize); 230 if (!y) SETERRQ(PETSC_COMM_SELF,1,"Interpolated y is not allocated"); 231 232 ierr = VecGetArray(X,&x_data);CHKERRQ(ierr); 233 N_VSetArrayPointer((realtype*)x_data,y); 234 ierr = CVodeGetDky(cvode->mem,t,0,y);CHKERRQ(ierr); 235 ierr = VecRestoreArray(X,&x_data);CHKERRQ(ierr); 236 PetscFunctionReturn(0); 237 } 238 239 #undef __FUNCT__ 240 #define __FUNCT__ "TSReset_Sundials" 241 PetscErrorCode TSReset_Sundials(TS ts) 242 { 243 TS_Sundials *cvode = (TS_Sundials*)ts->data; 244 PetscErrorCode ierr; 245 246 PetscFunctionBegin; 247 ierr = VecDestroy(&cvode->update);CHKERRQ(ierr); 248 ierr = VecDestroy(&cvode->ydot);CHKERRQ(ierr); 249 ierr = VecDestroy(&cvode->w1);CHKERRQ(ierr); 250 ierr = VecDestroy(&cvode->w2);CHKERRQ(ierr); 251 if (cvode->mem) CVodeFree(&cvode->mem); 252 PetscFunctionReturn(0); 253 } 254 255 #undef __FUNCT__ 256 #define __FUNCT__ "TSDestroy_Sundials" 257 PetscErrorCode TSDestroy_Sundials(TS ts) 258 { 259 TS_Sundials *cvode = (TS_Sundials*)ts->data; 260 PetscErrorCode ierr; 261 262 PetscFunctionBegin; 263 ierr = TSReset_Sundials(ts);CHKERRQ(ierr); 264 ierr = MPI_Comm_free(&(cvode->comm_sundials));CHKERRQ(ierr); 265 ierr = PetscFree(ts->data);CHKERRQ(ierr); 266 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetType_C","",PETSC_NULL);CHKERRQ(ierr); 267 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetMaxl_C","",PETSC_NULL);CHKERRQ(ierr); 268 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetLinearTolerance_C","",PETSC_NULL);CHKERRQ(ierr); 269 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetGramSchmidtType_C","",PETSC_NULL);CHKERRQ(ierr); 270 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetTolerance_C","",PETSC_NULL);CHKERRQ(ierr); 271 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetMinTimeStep_C","",PETSC_NULL);CHKERRQ(ierr); 272 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetMaxTimeStep_C","",PETSC_NULL);CHKERRQ(ierr); 273 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsGetPC_C","",PETSC_NULL);CHKERRQ(ierr); 274 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsGetIterations_C","",PETSC_NULL);CHKERRQ(ierr); 275 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsMonitorInternalSteps_C","",PETSC_NULL);CHKERRQ(ierr); 276 PetscFunctionReturn(0); 277 } 278 279 #undef __FUNCT__ 280 #define __FUNCT__ "TSSetUp_Sundials" 281 PetscErrorCode TSSetUp_Sundials(TS ts) 282 { 283 TS_Sundials *cvode = (TS_Sundials*)ts->data; 284 PetscErrorCode ierr; 285 PetscInt glosize,locsize,i,flag; 286 PetscScalar *y_data,*parray; 287 void *mem; 288 PC pc; 289 PCType pctype; 290 PetscBool pcnone; 291 292 PetscFunctionBegin; 293 /* get the vector size */ 294 ierr = VecGetSize(ts->vec_sol,&glosize);CHKERRQ(ierr); 295 ierr = VecGetLocalSize(ts->vec_sol,&locsize);CHKERRQ(ierr); 296 297 /* allocate the memory for N_Vec y */ 298 cvode->y = N_VNew_Parallel(cvode->comm_sundials,locsize,glosize); 299 if (!cvode->y) SETERRQ(PETSC_COMM_SELF,1,"cvode->y is not allocated"); 300 301 /* initialize N_Vec y: copy ts->vec_sol to cvode->y */ 302 ierr = VecGetArray(ts->vec_sol,&parray);CHKERRQ(ierr); 303 y_data = (PetscScalar*) N_VGetArrayPointer(cvode->y); 304 for (i = 0; i < locsize; i++) y_data[i] = parray[i]; 305 ierr = VecRestoreArray(ts->vec_sol,PETSC_NULL);CHKERRQ(ierr); 306 307 ierr = VecDuplicate(ts->vec_sol,&cvode->update);CHKERRQ(ierr); 308 ierr = VecDuplicate(ts->vec_sol,&cvode->ydot);CHKERRQ(ierr); 309 ierr = PetscLogObjectParent(ts,cvode->update);CHKERRQ(ierr); 310 ierr = PetscLogObjectParent(ts,cvode->ydot);CHKERRQ(ierr); 311 312 /* 313 Create work vectors for the TSPSolve_Sundials() routine. Note these are 314 allocated with zero space arrays because the actual array space is provided 315 by Sundials and set using VecPlaceArray(). 316 */ 317 ierr = VecCreateMPIWithArray(((PetscObject)ts)->comm,1,locsize,PETSC_DECIDE,0,&cvode->w1);CHKERRQ(ierr); 318 ierr = VecCreateMPIWithArray(((PetscObject)ts)->comm,1,locsize,PETSC_DECIDE,0,&cvode->w2);CHKERRQ(ierr); 319 ierr = PetscLogObjectParent(ts,cvode->w1);CHKERRQ(ierr); 320 ierr = PetscLogObjectParent(ts,cvode->w2);CHKERRQ(ierr); 321 322 /* Call CVodeCreate to create the solver memory and the use of a Newton iteration */ 323 mem = CVodeCreate(cvode->cvode_type, CV_NEWTON); 324 if (!mem) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"CVodeCreate() fails"); 325 cvode->mem = mem; 326 327 /* Set the pointer to user-defined data */ 328 flag = CVodeSetUserData(mem, ts); 329 if (flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSetUserData() fails"); 330 331 /* Sundials may choose to use a smaller initial step, but will never use a larger step. */ 332 flag = CVodeSetInitStep(mem,(realtype)ts->time_step); 333 if (flag) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_LIB,"CVodeSetInitStep() failed"); 334 if (cvode->mindt > 0) { 335 flag = CVodeSetMinStep(mem,(realtype)cvode->mindt); 336 if (flag) { 337 if (flag == CV_MEM_NULL) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_LIB,"CVodeSetMinStep() failed, cvode_mem pointer is NULL"); 338 else if (flag == CV_ILL_INPUT) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_LIB,"CVodeSetMinStep() failed, hmin is nonpositive or it exceeds the maximum allowable step size"); 339 else SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_LIB,"CVodeSetMinStep() failed"); 340 } 341 } 342 if (cvode->maxdt > 0) { 343 flag = CVodeSetMaxStep(mem,(realtype)cvode->maxdt); 344 if (flag) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_LIB,"CVodeSetMaxStep() failed"); 345 } 346 347 /* Call CVodeInit to initialize the integrator memory and specify the 348 * user's right hand side function in u'=f(t,u), the inital time T0, and 349 * the initial dependent variable vector cvode->y */ 350 flag = CVodeInit(mem,TSFunction_Sundials,ts->ptime,cvode->y); 351 if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeInit() fails, flag %d",flag); 352 353 /* specifies scalar relative and absolute tolerances */ 354 flag = CVodeSStolerances(mem,cvode->reltol,cvode->abstol); 355 if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSStolerances() fails, flag %d",flag); 356 357 /* Specify max num of steps to be taken by cvode in its attempt to reach the next output time */ 358 flag = CVodeSetMaxNumSteps(mem,ts->max_steps); 359 360 /* call CVSpgmr to use GMRES as the linear solver. */ 361 /* setup the ode integrator with the given preconditioner */ 362 ierr = TSSundialsGetPC(ts,&pc);CHKERRQ(ierr); 363 ierr = PCGetType(pc,&pctype);CHKERRQ(ierr); 364 ierr = PetscObjectTypeCompare((PetscObject)pc,PCNONE,&pcnone);CHKERRQ(ierr); 365 if (pcnone) { 366 flag = CVSpgmr(mem,PREC_NONE,0); 367 if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmr() fails, flag %d",flag); 368 } else { 369 flag = CVSpgmr(mem,PREC_LEFT,cvode->maxl); 370 if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmr() fails, flag %d",flag); 371 372 /* Set preconditioner and solve routines Precond and PSolve, 373 and the pointer to the user-defined block data */ 374 flag = CVSpilsSetPreconditioner(mem,TSPrecond_Sundials,TSPSolve_Sundials); 375 if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpilsSetPreconditioner() fails, flag %d", flag); 376 } 377 378 flag = CVSpilsSetGSType(mem, MODIFIED_GS); 379 if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmrSetGSType() fails, flag %d",flag); 380 PetscFunctionReturn(0); 381 } 382 383 /* type of CVODE linear multistep method */ 384 const char *const TSSundialsLmmTypes[] = {"","ADAMS","BDF","TSSundialsLmmType","SUNDIALS_",0}; 385 /* type of G-S orthogonalization used by CVODE linear solver */ 386 const char *const TSSundialsGramSchmidtTypes[] = {"","MODIFIED","CLASSICAL","TSSundialsGramSchmidtType","SUNDIALS_",0}; 387 388 #undef __FUNCT__ 389 #define __FUNCT__ "TSSetFromOptions_Sundials" 390 PetscErrorCode TSSetFromOptions_Sundials(TS ts) 391 { 392 TS_Sundials *cvode = (TS_Sundials*)ts->data; 393 PetscErrorCode ierr; 394 int indx; 395 PetscBool flag; 396 PC pc; 397 398 PetscFunctionBegin; 399 ierr = PetscOptionsHead("SUNDIALS ODE solver options");CHKERRQ(ierr); 400 ierr = PetscOptionsEList("-ts_sundials_type","Scheme","TSSundialsSetType",TSSundialsLmmTypes,3,TSSundialsLmmTypes[cvode->cvode_type],&indx,&flag);CHKERRQ(ierr); 401 if (flag) { 402 ierr = TSSundialsSetType(ts,(TSSundialsLmmType)indx);CHKERRQ(ierr); 403 } 404 ierr = PetscOptionsEList("-ts_sundials_gramschmidt_type","Type of orthogonalization","TSSundialsSetGramSchmidtType",TSSundialsGramSchmidtTypes,3,TSSundialsGramSchmidtTypes[cvode->gtype],&indx,&flag);CHKERRQ(ierr); 405 if (flag) { 406 ierr = TSSundialsSetGramSchmidtType(ts,(TSSundialsGramSchmidtType)indx);CHKERRQ(ierr); 407 } 408 ierr = PetscOptionsReal("-ts_sundials_atol","Absolute tolerance for convergence","TSSundialsSetTolerance",cvode->abstol,&cvode->abstol,PETSC_NULL);CHKERRQ(ierr); 409 ierr = PetscOptionsReal("-ts_sundials_rtol","Relative tolerance for convergence","TSSundialsSetTolerance",cvode->reltol,&cvode->reltol,PETSC_NULL);CHKERRQ(ierr); 410 ierr = PetscOptionsReal("-ts_sundials_mindt","Minimum step size","TSSundialsSetMinTimeStep",cvode->mindt,&cvode->mindt,PETSC_NULL);CHKERRQ(ierr); 411 ierr = PetscOptionsReal("-ts_sundials_maxdt","Maximum step size","TSSundialsSetMaxTimeStep",cvode->maxdt,&cvode->maxdt,PETSC_NULL);CHKERRQ(ierr); 412 ierr = PetscOptionsReal("-ts_sundials_linear_tolerance","Convergence tolerance for linear solve","TSSundialsSetLinearTolerance",cvode->linear_tol,&cvode->linear_tol,&flag);CHKERRQ(ierr); 413 ierr = PetscOptionsInt("-ts_sundials_maxl","Max dimension of the Krylov subspace","TSSundialsSetMaxl",cvode->maxl,&cvode->maxl,&flag);CHKERRQ(ierr); 414 ierr = PetscOptionsBool("-ts_sundials_monitor_steps","Monitor SUNDIALS internel steps","TSSundialsMonitorInternalSteps",cvode->monitorstep,&cvode->monitorstep,PETSC_NULL);CHKERRQ(ierr); 415 ierr = PetscOptionsTail();CHKERRQ(ierr); 416 ierr = TSSundialsGetPC(ts,&pc);CHKERRQ(ierr); 417 ierr = PCSetFromOptions(pc);CHKERRQ(ierr); 418 PetscFunctionReturn(0); 419 } 420 421 #undef __FUNCT__ 422 #define __FUNCT__ "TSView_Sundials" 423 PetscErrorCode TSView_Sundials(TS ts,PetscViewer viewer) 424 { 425 TS_Sundials *cvode = (TS_Sundials*)ts->data; 426 PetscErrorCode ierr; 427 char *type; 428 char atype[] = "Adams"; 429 char btype[] = "BDF: backward differentiation formula"; 430 PetscBool iascii,isstring; 431 long int nsteps,its,nfevals,nlinsetups,nfails,itmp; 432 PetscInt qlast,qcur; 433 PetscReal hinused,hlast,hcur,tcur,tolsfac; 434 PC pc; 435 436 PetscFunctionBegin; 437 if (cvode->cvode_type == SUNDIALS_ADAMS) type = atype; 438 else type = btype; 439 440 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 441 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 442 if (iascii) { 443 ierr = PetscViewerASCIIPrintf(viewer,"Sundials integrater does not use SNES!\n");CHKERRQ(ierr); 444 ierr = PetscViewerASCIIPrintf(viewer,"Sundials integrater type %s\n",type);CHKERRQ(ierr); 445 ierr = PetscViewerASCIIPrintf(viewer,"Sundials abs tol %g rel tol %g\n",cvode->abstol,cvode->reltol);CHKERRQ(ierr); 446 ierr = PetscViewerASCIIPrintf(viewer,"Sundials linear solver tolerance factor %g\n",cvode->linear_tol);CHKERRQ(ierr); 447 ierr = PetscViewerASCIIPrintf(viewer,"Sundials max dimension of Krylov subspace %D\n",cvode->maxl);CHKERRQ(ierr); 448 if (cvode->gtype == SUNDIALS_MODIFIED_GS) { 449 ierr = PetscViewerASCIIPrintf(viewer,"Sundials using modified Gram-Schmidt for orthogonalization in GMRES\n");CHKERRQ(ierr); 450 } else { 451 ierr = PetscViewerASCIIPrintf(viewer,"Sundials using unmodified (classical) Gram-Schmidt for orthogonalization in GMRES\n");CHKERRQ(ierr); 452 } 453 if (cvode->mindt > 0) {ierr = PetscViewerASCIIPrintf(viewer,"Sundials minimum time step %g\n",cvode->mindt);CHKERRQ(ierr);} 454 if (cvode->maxdt > 0) {ierr = PetscViewerASCIIPrintf(viewer,"Sundials maximum time step %g\n",cvode->maxdt);CHKERRQ(ierr);} 455 456 /* Outputs from CVODE, CVSPILS */ 457 ierr = CVodeGetTolScaleFactor(cvode->mem,&tolsfac);CHKERRQ(ierr); 458 ierr = PetscViewerASCIIPrintf(viewer,"Sundials suggested factor for tolerance scaling %g\n",tolsfac);CHKERRQ(ierr); 459 ierr = CVodeGetIntegratorStats(cvode->mem,&nsteps,&nfevals, 460 &nlinsetups,&nfails,&qlast,&qcur, 461 &hinused,&hlast,&hcur,&tcur);CHKERRQ(ierr); 462 ierr = PetscViewerASCIIPrintf(viewer,"Sundials cumulative number of internal steps %D\n",nsteps);CHKERRQ(ierr); 463 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to rhs function %D\n",nfevals);CHKERRQ(ierr); 464 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to linear solver setup function %D\n",nlinsetups);CHKERRQ(ierr); 465 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of error test failures %D\n",nfails);CHKERRQ(ierr); 466 467 ierr = CVodeGetNonlinSolvStats(cvode->mem,&its,&nfails);CHKERRQ(ierr); 468 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear solver iterations %D\n",its);CHKERRQ(ierr); 469 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear convergence failure %D\n",nfails);CHKERRQ(ierr); 470 471 ierr = CVSpilsGetNumLinIters(cvode->mem, &its);CHKERRQ(ierr); /* its = no. of calls to TSPrecond_Sundials() */ 472 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of linear iterations %D\n",its);CHKERRQ(ierr); 473 ierr = CVSpilsGetNumConvFails(cvode->mem,&itmp);CHKERRQ(ierr); 474 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of linear convergence failures %D\n",itmp);CHKERRQ(ierr); 475 476 ierr = TSSundialsGetPC(ts,&pc);CHKERRQ(ierr); 477 ierr = PCView(pc,viewer);CHKERRQ(ierr); 478 ierr = CVSpilsGetNumPrecEvals(cvode->mem,&itmp);CHKERRQ(ierr); 479 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner evaluations %D\n",itmp);CHKERRQ(ierr); 480 ierr = CVSpilsGetNumPrecSolves(cvode->mem,&itmp);CHKERRQ(ierr); 481 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner solves %D\n",itmp);CHKERRQ(ierr); 482 483 ierr = CVSpilsGetNumJtimesEvals(cvode->mem,&itmp);CHKERRQ(ierr); 484 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of Jacobian-vector product evaluations %D\n",itmp);CHKERRQ(ierr); 485 ierr = CVSpilsGetNumRhsEvals(cvode->mem,&itmp);CHKERRQ(ierr); 486 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of rhs calls for finite diff. Jacobian-vector evals %D\n",itmp);CHKERRQ(ierr); 487 } else if (isstring) { 488 ierr = PetscViewerStringSPrintf(viewer,"Sundials type %s",type);CHKERRQ(ierr); 489 } 490 PetscFunctionReturn(0); 491 } 492 493 494 /* --------------------------------------------------------------------------*/ 495 EXTERN_C_BEGIN 496 #undef __FUNCT__ 497 #define __FUNCT__ "TSSundialsSetType_Sundials" 498 PetscErrorCode TSSundialsSetType_Sundials(TS ts,TSSundialsLmmType type) 499 { 500 TS_Sundials *cvode = (TS_Sundials*)ts->data; 501 502 PetscFunctionBegin; 503 cvode->cvode_type = type; 504 PetscFunctionReturn(0); 505 } 506 EXTERN_C_END 507 508 EXTERN_C_BEGIN 509 #undef __FUNCT__ 510 #define __FUNCT__ "TSSundialsSetMaxl_Sundials" 511 PetscErrorCode TSSundialsSetMaxl_Sundials(TS ts,PetscInt maxl) 512 { 513 TS_Sundials *cvode = (TS_Sundials*)ts->data; 514 515 PetscFunctionBegin; 516 cvode->maxl = maxl; 517 PetscFunctionReturn(0); 518 } 519 EXTERN_C_END 520 521 EXTERN_C_BEGIN 522 #undef __FUNCT__ 523 #define __FUNCT__ "TSSundialsSetLinearTolerance_Sundials" 524 PetscErrorCode TSSundialsSetLinearTolerance_Sundials(TS ts,double tol) 525 { 526 TS_Sundials *cvode = (TS_Sundials*)ts->data; 527 528 PetscFunctionBegin; 529 cvode->linear_tol = tol; 530 PetscFunctionReturn(0); 531 } 532 EXTERN_C_END 533 534 EXTERN_C_BEGIN 535 #undef __FUNCT__ 536 #define __FUNCT__ "TSSundialsSetGramSchmidtType_Sundials" 537 PetscErrorCode TSSundialsSetGramSchmidtType_Sundials(TS ts,TSSundialsGramSchmidtType type) 538 { 539 TS_Sundials *cvode = (TS_Sundials*)ts->data; 540 541 PetscFunctionBegin; 542 cvode->gtype = type; 543 PetscFunctionReturn(0); 544 } 545 EXTERN_C_END 546 547 EXTERN_C_BEGIN 548 #undef __FUNCT__ 549 #define __FUNCT__ "TSSundialsSetTolerance_Sundials" 550 PetscErrorCode TSSundialsSetTolerance_Sundials(TS ts,double aabs,double rel) 551 { 552 TS_Sundials *cvode = (TS_Sundials*)ts->data; 553 554 PetscFunctionBegin; 555 if (aabs != PETSC_DECIDE) cvode->abstol = aabs; 556 if (rel != PETSC_DECIDE) cvode->reltol = rel; 557 PetscFunctionReturn(0); 558 } 559 EXTERN_C_END 560 561 EXTERN_C_BEGIN 562 #undef __FUNCT__ 563 #define __FUNCT__ "TSSundialsSetMinTimeStep_Sundials" 564 PetscErrorCode TSSundialsSetMinTimeStep_Sundials(TS ts,PetscReal mindt) 565 { 566 TS_Sundials *cvode = (TS_Sundials*)ts->data; 567 568 PetscFunctionBegin; 569 cvode->mindt = mindt; 570 PetscFunctionReturn(0); 571 } 572 EXTERN_C_END 573 574 EXTERN_C_BEGIN 575 #undef __FUNCT__ 576 #define __FUNCT__ "TSSundialsSetMaxTimeStep_Sundials" 577 PetscErrorCode TSSundialsSetMaxTimeStep_Sundials(TS ts,PetscReal maxdt) 578 { 579 TS_Sundials *cvode = (TS_Sundials*)ts->data; 580 581 PetscFunctionBegin; 582 cvode->maxdt = maxdt; 583 PetscFunctionReturn(0); 584 } 585 EXTERN_C_END 586 587 EXTERN_C_BEGIN 588 #undef __FUNCT__ 589 #define __FUNCT__ "TSSundialsGetPC_Sundials" 590 PetscErrorCode TSSundialsGetPC_Sundials(TS ts,PC *pc) 591 { 592 SNES snes; 593 KSP ksp; 594 PetscErrorCode ierr; 595 596 PetscFunctionBegin; 597 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 598 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 599 ierr = KSPGetPC(ksp,pc);CHKERRQ(ierr); 600 PetscFunctionReturn(0); 601 } 602 EXTERN_C_END 603 604 EXTERN_C_BEGIN 605 #undef __FUNCT__ 606 #define __FUNCT__ "TSSundialsGetIterations_Sundials" 607 PetscErrorCode TSSundialsGetIterations_Sundials(TS ts,int *nonlin,int *lin) 608 { 609 PetscFunctionBegin; 610 if (nonlin) *nonlin = ts->snes_its; 611 if (lin) *lin = ts->ksp_its; 612 PetscFunctionReturn(0); 613 } 614 EXTERN_C_END 615 616 EXTERN_C_BEGIN 617 #undef __FUNCT__ 618 #define __FUNCT__ "TSSundialsMonitorInternalSteps_Sundials" 619 PetscErrorCode TSSundialsMonitorInternalSteps_Sundials(TS ts,PetscBool s) 620 { 621 TS_Sundials *cvode = (TS_Sundials*)ts->data; 622 623 PetscFunctionBegin; 624 cvode->monitorstep = s; 625 PetscFunctionReturn(0); 626 } 627 EXTERN_C_END 628 /* -------------------------------------------------------------------------------------------*/ 629 630 #undef __FUNCT__ 631 #define __FUNCT__ "TSSundialsGetIterations" 632 /*@C 633 TSSundialsGetIterations - Gets the number of nonlinear and linear iterations used so far by Sundials. 634 635 Not Collective 636 637 Input parameters: 638 . ts - the time-step context 639 640 Output Parameters: 641 + nonlin - number of nonlinear iterations 642 - lin - number of linear iterations 643 644 Level: advanced 645 646 Notes: 647 These return the number since the creation of the TS object 648 649 .keywords: non-linear iterations, linear iterations 650 651 .seealso: TSSundialsSetType(), TSSundialsSetMaxl(), 652 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 653 TSSundialsGetIterations(), TSSundialsSetType(), 654 TSSundialsSetLinearTolerance(), TSSundialsGetPC(), TSSetExactFinalTime() 655 656 @*/ 657 PetscErrorCode TSSundialsGetIterations(TS ts,int *nonlin,int *lin) 658 { 659 PetscErrorCode ierr; 660 661 PetscFunctionBegin; 662 ierr = PetscUseMethod(ts,"TSSundialsGetIterations_C",(TS,int*,int*),(ts,nonlin,lin));CHKERRQ(ierr); 663 PetscFunctionReturn(0); 664 } 665 666 #undef __FUNCT__ 667 #define __FUNCT__ "TSSundialsSetType" 668 /*@ 669 TSSundialsSetType - Sets the method that Sundials will use for integration. 670 671 Logically Collective on TS 672 673 Input parameters: 674 + ts - the time-step context 675 - type - one of SUNDIALS_ADAMS or SUNDIALS_BDF 676 677 Level: intermediate 678 679 .keywords: Adams, backward differentiation formula 680 681 .seealso: TSSundialsGetIterations(), TSSundialsSetMaxl(), 682 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 683 TSSundialsGetIterations(), TSSundialsSetType(), 684 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(), 685 TSSetExactFinalTime() 686 @*/ 687 PetscErrorCode TSSundialsSetType(TS ts,TSSundialsLmmType type) 688 { 689 PetscErrorCode ierr; 690 691 PetscFunctionBegin; 692 ierr = PetscTryMethod(ts,"TSSundialsSetType_C",(TS,TSSundialsLmmType),(ts,type));CHKERRQ(ierr); 693 PetscFunctionReturn(0); 694 } 695 696 #undef __FUNCT__ 697 #define __FUNCT__ "TSSundialsSetMaxl" 698 /*@ 699 TSSundialsSetMaxl - Sets the dimension of the Krylov space used by 700 GMRES in the linear solver in SUNDIALS. SUNDIALS DOES NOT use restarted GMRES so 701 this is the maximum number of GMRES steps that will be used. 702 703 Logically Collective on TS 704 705 Input parameters: 706 + ts - the time-step context 707 - maxl - number of direction vectors (the dimension of Krylov subspace). 708 709 Level: advanced 710 711 .keywords: GMRES 712 713 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), 714 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 715 TSSundialsGetIterations(), TSSundialsSetType(), 716 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(), 717 TSSetExactFinalTime() 718 719 @*/ 720 PetscErrorCode TSSundialsSetMaxl(TS ts,PetscInt maxl) 721 { 722 PetscErrorCode ierr; 723 724 PetscFunctionBegin; 725 PetscValidLogicalCollectiveInt(ts,maxl,2); 726 ierr = PetscTryMethod(ts,"TSSundialsSetMaxl_C",(TS,PetscInt),(ts,maxl));CHKERRQ(ierr); 727 PetscFunctionReturn(0); 728 } 729 730 #undef __FUNCT__ 731 #define __FUNCT__ "TSSundialsSetLinearTolerance" 732 /*@ 733 TSSundialsSetLinearTolerance - Sets the tolerance used to solve the linear 734 system by SUNDIALS. 735 736 Logically Collective on TS 737 738 Input parameters: 739 + ts - the time-step context 740 - tol - the factor by which the tolerance on the nonlinear solver is 741 multiplied to get the tolerance on the linear solver, .05 by default. 742 743 Level: advanced 744 745 .keywords: GMRES, linear convergence tolerance, SUNDIALS 746 747 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(), 748 TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 749 TSSundialsGetIterations(), TSSundialsSetType(), 750 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(), 751 TSSetExactFinalTime() 752 753 @*/ 754 PetscErrorCode TSSundialsSetLinearTolerance(TS ts,double tol) 755 { 756 PetscErrorCode ierr; 757 758 PetscFunctionBegin; 759 PetscValidLogicalCollectiveReal(ts,tol,2); 760 ierr = PetscTryMethod(ts,"TSSundialsSetLinearTolerance_C",(TS,double),(ts,tol));CHKERRQ(ierr); 761 PetscFunctionReturn(0); 762 } 763 764 #undef __FUNCT__ 765 #define __FUNCT__ "TSSundialsSetGramSchmidtType" 766 /*@ 767 TSSundialsSetGramSchmidtType - Sets type of orthogonalization used 768 in GMRES method by SUNDIALS linear solver. 769 770 Logically Collective on TS 771 772 Input parameters: 773 + ts - the time-step context 774 - type - either SUNDIALS_MODIFIED_GS or SUNDIALS_CLASSICAL_GS 775 776 Level: advanced 777 778 .keywords: Sundials, orthogonalization 779 780 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(), 781 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), 782 TSSundialsGetIterations(), TSSundialsSetType(), 783 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(), 784 TSSetExactFinalTime() 785 786 @*/ 787 PetscErrorCode TSSundialsSetGramSchmidtType(TS ts,TSSundialsGramSchmidtType type) 788 { 789 PetscErrorCode ierr; 790 791 PetscFunctionBegin; 792 ierr = PetscTryMethod(ts,"TSSundialsSetGramSchmidtType_C",(TS,TSSundialsGramSchmidtType),(ts,type));CHKERRQ(ierr); 793 PetscFunctionReturn(0); 794 } 795 796 #undef __FUNCT__ 797 #define __FUNCT__ "TSSundialsSetTolerance" 798 /*@ 799 TSSundialsSetTolerance - Sets the absolute and relative tolerance used by 800 Sundials for error control. 801 802 Logically Collective on TS 803 804 Input parameters: 805 + ts - the time-step context 806 . aabs - the absolute tolerance 807 - rel - the relative tolerance 808 809 See the Cvode/Sundials users manual for exact details on these parameters. Essentially 810 these regulate the size of the error for a SINGLE timestep. 811 812 Level: intermediate 813 814 .keywords: Sundials, tolerance 815 816 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESMaxl(), 817 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), 818 TSSundialsGetIterations(), TSSundialsSetType(), 819 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(), 820 TSSetExactFinalTime() 821 822 @*/ 823 PetscErrorCode TSSundialsSetTolerance(TS ts,double aabs,double rel) 824 { 825 PetscErrorCode ierr; 826 827 PetscFunctionBegin; 828 ierr = PetscTryMethod(ts,"TSSundialsSetTolerance_C",(TS,double,double),(ts,aabs,rel));CHKERRQ(ierr); 829 PetscFunctionReturn(0); 830 } 831 832 #undef __FUNCT__ 833 #define __FUNCT__ "TSSundialsGetPC" 834 /*@ 835 TSSundialsGetPC - Extract the PC context from a time-step context for Sundials. 836 837 Input Parameter: 838 . ts - the time-step context 839 840 Output Parameter: 841 . pc - the preconditioner context 842 843 Level: advanced 844 845 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(), 846 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 847 TSSundialsGetIterations(), TSSundialsSetType(), 848 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance() 849 @*/ 850 PetscErrorCode TSSundialsGetPC(TS ts,PC *pc) 851 { 852 PetscErrorCode ierr; 853 854 PetscFunctionBegin; 855 ierr = PetscUseMethod(ts,"TSSundialsGetPC_C",(TS,PC*),(ts,pc));CHKERRQ(ierr); 856 PetscFunctionReturn(0); 857 } 858 859 #undef __FUNCT__ 860 #define __FUNCT__ "TSSundialsSetMinTimeStep" 861 /*@ 862 TSSundialsSetMinTimeStep - Smallest time step to be chosen by the adaptive controller. 863 864 Input Parameter: 865 + ts - the time-step context 866 - mindt - lowest time step if positive, negative to deactivate 867 868 Note: 869 Sundials will error if it is not possible to keep the estimated truncation error below 870 the tolerance set with TSSundialsSetTolerance() without going below this step size. 871 872 Level: beginner 873 874 .seealso: TSSundialsSetType(), TSSundialsSetTolerance(), 875 @*/ 876 PetscErrorCode TSSundialsSetMinTimeStep(TS ts,PetscReal mindt) 877 { 878 PetscErrorCode ierr; 879 880 PetscFunctionBegin; 881 ierr = PetscTryMethod(ts,"TSSundialsSetMinTimeStep_C",(TS,PetscReal),(ts,mindt));CHKERRQ(ierr); 882 PetscFunctionReturn(0); 883 } 884 885 #undef __FUNCT__ 886 #define __FUNCT__ "TSSundialsSetMaxTimeStep" 887 /*@ 888 TSSundialsSetMaxTimeStep - Largest time step to be chosen by the adaptive controller. 889 890 Input Parameter: 891 + ts - the time-step context 892 - maxdt - lowest time step if positive, negative to deactivate 893 894 Level: beginner 895 896 .seealso: TSSundialsSetType(), TSSundialsSetTolerance(), 897 @*/ 898 PetscErrorCode TSSundialsSetMaxTimeStep(TS ts,PetscReal maxdt) 899 { 900 PetscErrorCode ierr; 901 902 PetscFunctionBegin; 903 ierr = PetscTryMethod(ts,"TSSundialsSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));CHKERRQ(ierr); 904 PetscFunctionReturn(0); 905 } 906 907 #undef __FUNCT__ 908 #define __FUNCT__ "TSSundialsMonitorInternalSteps" 909 /*@ 910 TSSundialsMonitorInternalSteps - Monitor Sundials internal steps (Defaults to false). 911 912 Input Parameter: 913 + ts - the time-step context 914 - ft - PETSC_TRUE if monitor, else PETSC_FALSE 915 916 Level: beginner 917 918 .seealso:TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(), 919 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 920 TSSundialsGetIterations(), TSSundialsSetType(), 921 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC() 922 @*/ 923 PetscErrorCode TSSundialsMonitorInternalSteps(TS ts,PetscBool ft) 924 { 925 PetscErrorCode ierr; 926 927 PetscFunctionBegin; 928 ierr = PetscTryMethod(ts,"TSSundialsMonitorInternalSteps_C",(TS,PetscBool),(ts,ft));CHKERRQ(ierr); 929 PetscFunctionReturn(0); 930 } 931 /* -------------------------------------------------------------------------------------------*/ 932 /*MC 933 TSSUNDIALS - ODE solver using the LLNL CVODE/SUNDIALS package (now called SUNDIALS) 934 935 Options Database: 936 + -ts_sundials_type <bdf,adams> 937 . -ts_sundials_gramschmidt_type <modified, classical> - type of orthogonalization inside GMRES 938 . -ts_sundials_atol <tol> - Absolute tolerance for convergence 939 . -ts_sundials_rtol <tol> - Relative tolerance for convergence 940 . -ts_sundials_linear_tolerance <tol> 941 . -ts_sundials_maxl <maxl> - Max dimension of the Krylov subspace 942 - -ts_sundials_monitor_steps - Monitor SUNDIALS internel steps 943 944 945 Notes: This uses its own nonlinear solver and Krylov method so PETSc SNES and KSP options do not apply 946 only PETSc PC options 947 948 Level: beginner 949 950 .seealso: TSCreate(), TS, TSSetType(), TSSundialsSetType(), TSSundialsSetMaxl(), TSSundialsSetLinearTolerance(), 951 TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), TSSundialsGetPC(), TSSundialsGetIterations(), TSSetExactFinalTime() 952 953 M*/ 954 EXTERN_C_BEGIN 955 #undef __FUNCT__ 956 #define __FUNCT__ "TSCreate_Sundials" 957 PetscErrorCode TSCreate_Sundials(TS ts) 958 { 959 TS_Sundials *cvode; 960 PetscErrorCode ierr; 961 PC pc; 962 963 PetscFunctionBegin; 964 ts->ops->reset = TSReset_Sundials; 965 ts->ops->destroy = TSDestroy_Sundials; 966 ts->ops->view = TSView_Sundials; 967 ts->ops->setup = TSSetUp_Sundials; 968 ts->ops->step = TSStep_Sundials; 969 ts->ops->interpolate = TSInterpolate_Sundials; 970 ts->ops->setfromoptions = TSSetFromOptions_Sundials; 971 972 ierr = PetscNewLog(ts,TS_Sundials,&cvode);CHKERRQ(ierr); 973 974 ts->data = (void*)cvode; 975 cvode->cvode_type = SUNDIALS_BDF; 976 cvode->gtype = SUNDIALS_CLASSICAL_GS; 977 cvode->maxl = 5; 978 cvode->linear_tol = .05; 979 cvode->monitorstep = PETSC_TRUE; 980 981 ierr = MPI_Comm_dup(((PetscObject)ts)->comm,&(cvode->comm_sundials));CHKERRQ(ierr); 982 983 cvode->mindt = -1.; 984 cvode->maxdt = -1.; 985 986 /* set tolerance for Sundials */ 987 cvode->reltol = 1e-6; 988 cvode->abstol = 1e-6; 989 990 /* set PCNONE as default pctype */ 991 ierr = TSSundialsGetPC_Sundials(ts,&pc);CHKERRQ(ierr); 992 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 993 994 if (ts->exact_final_time != TS_EXACTFINALTIME_STEPOVER) ts->exact_final_time = TS_EXACTFINALTIME_STEPOVER; 995 996 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetType_C","TSSundialsSetType_Sundials", 997 TSSundialsSetType_Sundials);CHKERRQ(ierr); 998 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetMaxl_C", 999 "TSSundialsSetMaxl_Sundials", 1000 TSSundialsSetMaxl_Sundials);CHKERRQ(ierr); 1001 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetLinearTolerance_C", 1002 "TSSundialsSetLinearTolerance_Sundials", 1003 TSSundialsSetLinearTolerance_Sundials);CHKERRQ(ierr); 1004 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetGramSchmidtType_C", 1005 "TSSundialsSetGramSchmidtType_Sundials", 1006 TSSundialsSetGramSchmidtType_Sundials);CHKERRQ(ierr); 1007 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetTolerance_C", 1008 "TSSundialsSetTolerance_Sundials", 1009 TSSundialsSetTolerance_Sundials);CHKERRQ(ierr); 1010 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetMinTimeStep_C", 1011 "TSSundialsSetMinTimeStep_Sundials", 1012 TSSundialsSetMinTimeStep_Sundials);CHKERRQ(ierr); 1013 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetMaxTimeStep_C", 1014 "TSSundialsSetMaxTimeStep_Sundials", 1015 TSSundialsSetMaxTimeStep_Sundials);CHKERRQ(ierr); 1016 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsGetPC_C", 1017 "TSSundialsGetPC_Sundials", 1018 TSSundialsGetPC_Sundials);CHKERRQ(ierr); 1019 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsGetIterations_C", 1020 "TSSundialsGetIterations_Sundials", 1021 TSSundialsGetIterations_Sundials);CHKERRQ(ierr); 1022 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsMonitorInternalSteps_C", 1023 "TSSundialsMonitorInternalSteps_Sundials", 1024 TSSundialsMonitorInternalSteps_Sundials);CHKERRQ(ierr); 1025 PetscFunctionReturn(0); 1026 } 1027 EXTERN_C_END 1028