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 <../src/ts/impls/implicit/sundials/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 PetscScalar *y_data; 27 28 PetscFunctionBegin; 29 ierr = TSGetIJacobian(ts,&J,&P,NULL,NULL);CHKERRQ(ierr); 30 y_data = (PetscScalar*) N_VGetArrayPointer(y); 31 ierr = VecPlaceArray(yy,y_data);CHKERRQ(ierr); 32 ierr = VecZeroEntries(yydot);CHKERRQ(ierr); /* The Jacobian is independent of Ydot for ODE which is all that CVode works for */ 33 /* compute the shifted Jacobian (1/gm)*I + Jrest */ 34 ierr = TSComputeIJacobian(ts,ts->ptime,yy,yydot,1/gm,J,P,PETSC_FALSE);CHKERRQ(ierr); 35 ierr = VecResetArray(yy);CHKERRQ(ierr); 36 ierr = MatScale(P,gm);CHKERRQ(ierr); /* turn into I-gm*Jrest, J is not used by Sundials */ 37 *jcurPtr = TRUE; 38 ierr = TSSundialsGetPC(ts,&pc);CHKERRQ(ierr); 39 ierr = PCSetOperators(pc,J,P);CHKERRQ(ierr); 40 PetscFunctionReturn(0); 41 } 42 43 /* 44 TSPSolve_Sundials - routine that we provide to Sundials that applies the preconditioner. 45 */ 46 #undef __FUNCT__ 47 #define __FUNCT__ "TSPSolve_Sundials" 48 PetscErrorCode TSPSolve_Sundials(realtype tn,N_Vector y,N_Vector fy,N_Vector r,N_Vector z, 49 realtype _gamma,realtype delta,int lr,void *P_data,N_Vector vtemp) 50 { 51 TS ts = (TS) P_data; 52 TS_Sundials *cvode = (TS_Sundials*)ts->data; 53 PC pc; 54 Vec rr = cvode->w1,zz = cvode->w2; 55 PetscErrorCode ierr; 56 PetscScalar *r_data,*z_data; 57 58 PetscFunctionBegin; 59 /* Make the PETSc work vectors rr and zz point to the arrays in the SUNDIALS vectors r and z respectively*/ 60 r_data = (PetscScalar*) N_VGetArrayPointer(r); 61 z_data = (PetscScalar*) N_VGetArrayPointer(z); 62 ierr = VecPlaceArray(rr,r_data);CHKERRQ(ierr); 63 ierr = VecPlaceArray(zz,z_data);CHKERRQ(ierr); 64 65 /* Solve the Px=r and put the result in zz */ 66 ierr = TSSundialsGetPC(ts,&pc);CHKERRQ(ierr); 67 ierr = PCApply(pc,rr,zz);CHKERRQ(ierr); 68 ierr = VecResetArray(rr);CHKERRQ(ierr); 69 ierr = VecResetArray(zz);CHKERRQ(ierr); 70 PetscFunctionReturn(0); 71 } 72 73 /* 74 TSFunction_Sundials - routine that we provide to Sundials that applies the right hand side. 75 */ 76 #undef __FUNCT__ 77 #define __FUNCT__ "TSFunction_Sundials" 78 int TSFunction_Sundials(realtype t,N_Vector y,N_Vector ydot,void *ctx) 79 { 80 TS ts = (TS) ctx; 81 DM dm; 82 DMTS tsdm; 83 TSIFunction ifunction; 84 MPI_Comm comm; 85 TS_Sundials *cvode = (TS_Sundials*)ts->data; 86 Vec yy = cvode->w1,yyd = cvode->w2,yydot = cvode->ydot; 87 PetscScalar *y_data,*ydot_data; 88 PetscErrorCode ierr; 89 90 PetscFunctionBegin; 91 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 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,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,NULL);CHKERRQ(ierr); 135 136 ierr = TSPreStep(ts);CHKERRQ(ierr); 137 138 /* We would like to call TSPreStep() when starting each step (including rejections), TSPreStage(), 139 * and TSPostStage() before each 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()",(double)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 = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",NULL);CHKERRQ(ierr); 267 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",NULL);CHKERRQ(ierr); 268 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",NULL);CHKERRQ(ierr); 269 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",NULL);CHKERRQ(ierr); 270 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",NULL);CHKERRQ(ierr); 271 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",NULL);CHKERRQ(ierr); 272 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",NULL);CHKERRQ(ierr); 273 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",NULL);CHKERRQ(ierr); 274 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",NULL);CHKERRQ(ierr); 275 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",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,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((PetscObject)ts,(PetscObject)cvode->update);CHKERRQ(ierr); 310 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)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(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,0,&cvode->w1);CHKERRQ(ierr); 318 ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,0,&cvode->w2);CHKERRQ(ierr); 319 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w1);CHKERRQ(ierr); 320 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)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(PetscObjectComm((PetscObject)ts),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(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed, cvode_mem pointer is NULL"); 338 else if (flag == CV_ILL_INPUT) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed, hmin is nonpositive or it exceeds the maximum allowable step size"); 339 else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed"); 340 } 341 } 342 if (cvode->maxdt > 0) { 343 flag = CVodeSetMaxStep(mem,(realtype)cvode->maxdt); 344 if (flag) SETERRQ(PetscObjectComm((PetscObject)ts),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,NULL);CHKERRQ(ierr); 409 ierr = PetscOptionsReal("-ts_sundials_rtol","Relative tolerance for convergence","TSSundialsSetTolerance",cvode->reltol,&cvode->reltol,NULL);CHKERRQ(ierr); 410 ierr = PetscOptionsReal("-ts_sundials_mindt","Minimum step size","TSSundialsSetMinTimeStep",cvode->mindt,&cvode->mindt,NULL);CHKERRQ(ierr); 411 ierr = PetscOptionsReal("-ts_sundials_maxdt","Maximum step size","TSSundialsSetMaxTimeStep",cvode->maxdt,&cvode->maxdt,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,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 #undef __FUNCT__ 496 #define __FUNCT__ "TSSundialsSetType_Sundials" 497 PetscErrorCode TSSundialsSetType_Sundials(TS ts,TSSundialsLmmType type) 498 { 499 TS_Sundials *cvode = (TS_Sundials*)ts->data; 500 501 PetscFunctionBegin; 502 cvode->cvode_type = type; 503 PetscFunctionReturn(0); 504 } 505 506 #undef __FUNCT__ 507 #define __FUNCT__ "TSSundialsSetMaxl_Sundials" 508 PetscErrorCode TSSundialsSetMaxl_Sundials(TS ts,PetscInt maxl) 509 { 510 TS_Sundials *cvode = (TS_Sundials*)ts->data; 511 512 PetscFunctionBegin; 513 cvode->maxl = maxl; 514 PetscFunctionReturn(0); 515 } 516 517 #undef __FUNCT__ 518 #define __FUNCT__ "TSSundialsSetLinearTolerance_Sundials" 519 PetscErrorCode TSSundialsSetLinearTolerance_Sundials(TS ts,double tol) 520 { 521 TS_Sundials *cvode = (TS_Sundials*)ts->data; 522 523 PetscFunctionBegin; 524 cvode->linear_tol = tol; 525 PetscFunctionReturn(0); 526 } 527 528 #undef __FUNCT__ 529 #define __FUNCT__ "TSSundialsSetGramSchmidtType_Sundials" 530 PetscErrorCode TSSundialsSetGramSchmidtType_Sundials(TS ts,TSSundialsGramSchmidtType type) 531 { 532 TS_Sundials *cvode = (TS_Sundials*)ts->data; 533 534 PetscFunctionBegin; 535 cvode->gtype = type; 536 PetscFunctionReturn(0); 537 } 538 539 #undef __FUNCT__ 540 #define __FUNCT__ "TSSundialsSetTolerance_Sundials" 541 PetscErrorCode TSSundialsSetTolerance_Sundials(TS ts,double aabs,double rel) 542 { 543 TS_Sundials *cvode = (TS_Sundials*)ts->data; 544 545 PetscFunctionBegin; 546 if (aabs != PETSC_DECIDE) cvode->abstol = aabs; 547 if (rel != PETSC_DECIDE) cvode->reltol = rel; 548 PetscFunctionReturn(0); 549 } 550 551 #undef __FUNCT__ 552 #define __FUNCT__ "TSSundialsSetMinTimeStep_Sundials" 553 PetscErrorCode TSSundialsSetMinTimeStep_Sundials(TS ts,PetscReal mindt) 554 { 555 TS_Sundials *cvode = (TS_Sundials*)ts->data; 556 557 PetscFunctionBegin; 558 cvode->mindt = mindt; 559 PetscFunctionReturn(0); 560 } 561 562 #undef __FUNCT__ 563 #define __FUNCT__ "TSSundialsSetMaxTimeStep_Sundials" 564 PetscErrorCode TSSundialsSetMaxTimeStep_Sundials(TS ts,PetscReal maxdt) 565 { 566 TS_Sundials *cvode = (TS_Sundials*)ts->data; 567 568 PetscFunctionBegin; 569 cvode->maxdt = maxdt; 570 PetscFunctionReturn(0); 571 } 572 #undef __FUNCT__ 573 #define __FUNCT__ "TSSundialsGetPC_Sundials" 574 PetscErrorCode TSSundialsGetPC_Sundials(TS ts,PC *pc) 575 { 576 SNES snes; 577 KSP ksp; 578 PetscErrorCode ierr; 579 580 PetscFunctionBegin; 581 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 582 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 583 ierr = KSPGetPC(ksp,pc);CHKERRQ(ierr); 584 PetscFunctionReturn(0); 585 } 586 587 #undef __FUNCT__ 588 #define __FUNCT__ "TSSundialsGetIterations_Sundials" 589 PetscErrorCode TSSundialsGetIterations_Sundials(TS ts,int *nonlin,int *lin) 590 { 591 PetscFunctionBegin; 592 if (nonlin) *nonlin = ts->snes_its; 593 if (lin) *lin = ts->ksp_its; 594 PetscFunctionReturn(0); 595 } 596 597 #undef __FUNCT__ 598 #define __FUNCT__ "TSSundialsMonitorInternalSteps_Sundials" 599 PetscErrorCode TSSundialsMonitorInternalSteps_Sundials(TS ts,PetscBool s) 600 { 601 TS_Sundials *cvode = (TS_Sundials*)ts->data; 602 603 PetscFunctionBegin; 604 cvode->monitorstep = s; 605 PetscFunctionReturn(0); 606 } 607 /* -------------------------------------------------------------------------------------------*/ 608 609 #undef __FUNCT__ 610 #define __FUNCT__ "TSSundialsGetIterations" 611 /*@C 612 TSSundialsGetIterations - Gets the number of nonlinear and linear iterations used so far by Sundials. 613 614 Not Collective 615 616 Input parameters: 617 . ts - the time-step context 618 619 Output Parameters: 620 + nonlin - number of nonlinear iterations 621 - lin - number of linear iterations 622 623 Level: advanced 624 625 Notes: 626 These return the number since the creation of the TS object 627 628 .keywords: non-linear iterations, linear iterations 629 630 .seealso: TSSundialsSetType(), TSSundialsSetMaxl(), 631 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 632 TSSundialsGetIterations(), TSSundialsSetType(), 633 TSSundialsSetLinearTolerance(), TSSundialsGetPC(), TSSetExactFinalTime() 634 635 @*/ 636 PetscErrorCode TSSundialsGetIterations(TS ts,int *nonlin,int *lin) 637 { 638 PetscErrorCode ierr; 639 640 PetscFunctionBegin; 641 ierr = PetscUseMethod(ts,"TSSundialsGetIterations_C",(TS,int*,int*),(ts,nonlin,lin));CHKERRQ(ierr); 642 PetscFunctionReturn(0); 643 } 644 645 #undef __FUNCT__ 646 #define __FUNCT__ "TSSundialsSetType" 647 /*@ 648 TSSundialsSetType - Sets the method that Sundials will use for integration. 649 650 Logically Collective on TS 651 652 Input parameters: 653 + ts - the time-step context 654 - type - one of SUNDIALS_ADAMS or SUNDIALS_BDF 655 656 Level: intermediate 657 658 .keywords: Adams, backward differentiation formula 659 660 .seealso: TSSundialsGetIterations(), TSSundialsSetMaxl(), 661 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 662 TSSundialsGetIterations(), TSSundialsSetType(), 663 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(), 664 TSSetExactFinalTime() 665 @*/ 666 PetscErrorCode TSSundialsSetType(TS ts,TSSundialsLmmType type) 667 { 668 PetscErrorCode ierr; 669 670 PetscFunctionBegin; 671 ierr = PetscTryMethod(ts,"TSSundialsSetType_C",(TS,TSSundialsLmmType),(ts,type));CHKERRQ(ierr); 672 PetscFunctionReturn(0); 673 } 674 675 #undef __FUNCT__ 676 #define __FUNCT__ "TSSundialsSetMaxl" 677 /*@ 678 TSSundialsSetMaxl - Sets the dimension of the Krylov space used by 679 GMRES in the linear solver in SUNDIALS. SUNDIALS DOES NOT use restarted GMRES so 680 this is the maximum number of GMRES steps that will be used. 681 682 Logically Collective on TS 683 684 Input parameters: 685 + ts - the time-step context 686 - maxl - number of direction vectors (the dimension of Krylov subspace). 687 688 Level: advanced 689 690 .keywords: GMRES 691 692 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), 693 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 694 TSSundialsGetIterations(), TSSundialsSetType(), 695 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(), 696 TSSetExactFinalTime() 697 698 @*/ 699 PetscErrorCode TSSundialsSetMaxl(TS ts,PetscInt maxl) 700 { 701 PetscErrorCode ierr; 702 703 PetscFunctionBegin; 704 PetscValidLogicalCollectiveInt(ts,maxl,2); 705 ierr = PetscTryMethod(ts,"TSSundialsSetMaxl_C",(TS,PetscInt),(ts,maxl));CHKERRQ(ierr); 706 PetscFunctionReturn(0); 707 } 708 709 #undef __FUNCT__ 710 #define __FUNCT__ "TSSundialsSetLinearTolerance" 711 /*@ 712 TSSundialsSetLinearTolerance - Sets the tolerance used to solve the linear 713 system by SUNDIALS. 714 715 Logically Collective on TS 716 717 Input parameters: 718 + ts - the time-step context 719 - tol - the factor by which the tolerance on the nonlinear solver is 720 multiplied to get the tolerance on the linear solver, .05 by default. 721 722 Level: advanced 723 724 .keywords: GMRES, linear convergence tolerance, SUNDIALS 725 726 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(), 727 TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 728 TSSundialsGetIterations(), TSSundialsSetType(), 729 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(), 730 TSSetExactFinalTime() 731 732 @*/ 733 PetscErrorCode TSSundialsSetLinearTolerance(TS ts,double tol) 734 { 735 PetscErrorCode ierr; 736 737 PetscFunctionBegin; 738 PetscValidLogicalCollectiveReal(ts,tol,2); 739 ierr = PetscTryMethod(ts,"TSSundialsSetLinearTolerance_C",(TS,double),(ts,tol));CHKERRQ(ierr); 740 PetscFunctionReturn(0); 741 } 742 743 #undef __FUNCT__ 744 #define __FUNCT__ "TSSundialsSetGramSchmidtType" 745 /*@ 746 TSSundialsSetGramSchmidtType - Sets type of orthogonalization used 747 in GMRES method by SUNDIALS linear solver. 748 749 Logically Collective on TS 750 751 Input parameters: 752 + ts - the time-step context 753 - type - either SUNDIALS_MODIFIED_GS or SUNDIALS_CLASSICAL_GS 754 755 Level: advanced 756 757 .keywords: Sundials, orthogonalization 758 759 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(), 760 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), 761 TSSundialsGetIterations(), TSSundialsSetType(), 762 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(), 763 TSSetExactFinalTime() 764 765 @*/ 766 PetscErrorCode TSSundialsSetGramSchmidtType(TS ts,TSSundialsGramSchmidtType type) 767 { 768 PetscErrorCode ierr; 769 770 PetscFunctionBegin; 771 ierr = PetscTryMethod(ts,"TSSundialsSetGramSchmidtType_C",(TS,TSSundialsGramSchmidtType),(ts,type));CHKERRQ(ierr); 772 PetscFunctionReturn(0); 773 } 774 775 #undef __FUNCT__ 776 #define __FUNCT__ "TSSundialsSetTolerance" 777 /*@ 778 TSSundialsSetTolerance - Sets the absolute and relative tolerance used by 779 Sundials for error control. 780 781 Logically Collective on TS 782 783 Input parameters: 784 + ts - the time-step context 785 . aabs - the absolute tolerance 786 - rel - the relative tolerance 787 788 See the Cvode/Sundials users manual for exact details on these parameters. Essentially 789 these regulate the size of the error for a SINGLE timestep. 790 791 Level: intermediate 792 793 .keywords: Sundials, tolerance 794 795 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESMaxl(), 796 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), 797 TSSundialsGetIterations(), TSSundialsSetType(), 798 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(), 799 TSSetExactFinalTime() 800 801 @*/ 802 PetscErrorCode TSSundialsSetTolerance(TS ts,double aabs,double rel) 803 { 804 PetscErrorCode ierr; 805 806 PetscFunctionBegin; 807 ierr = PetscTryMethod(ts,"TSSundialsSetTolerance_C",(TS,double,double),(ts,aabs,rel));CHKERRQ(ierr); 808 PetscFunctionReturn(0); 809 } 810 811 #undef __FUNCT__ 812 #define __FUNCT__ "TSSundialsGetPC" 813 /*@ 814 TSSundialsGetPC - Extract the PC context from a time-step context for Sundials. 815 816 Input Parameter: 817 . ts - the time-step context 818 819 Output Parameter: 820 . pc - the preconditioner context 821 822 Level: advanced 823 824 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(), 825 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 826 TSSundialsGetIterations(), TSSundialsSetType(), 827 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance() 828 @*/ 829 PetscErrorCode TSSundialsGetPC(TS ts,PC *pc) 830 { 831 PetscErrorCode ierr; 832 833 PetscFunctionBegin; 834 ierr = PetscUseMethod(ts,"TSSundialsGetPC_C",(TS,PC*),(ts,pc));CHKERRQ(ierr); 835 PetscFunctionReturn(0); 836 } 837 838 #undef __FUNCT__ 839 #define __FUNCT__ "TSSundialsSetMinTimeStep" 840 /*@ 841 TSSundialsSetMinTimeStep - Smallest time step to be chosen by the adaptive controller. 842 843 Input Parameter: 844 + ts - the time-step context 845 - mindt - lowest time step if positive, negative to deactivate 846 847 Note: 848 Sundials will error if it is not possible to keep the estimated truncation error below 849 the tolerance set with TSSundialsSetTolerance() without going below this step size. 850 851 Level: beginner 852 853 .seealso: TSSundialsSetType(), TSSundialsSetTolerance(), 854 @*/ 855 PetscErrorCode TSSundialsSetMinTimeStep(TS ts,PetscReal mindt) 856 { 857 PetscErrorCode ierr; 858 859 PetscFunctionBegin; 860 ierr = PetscTryMethod(ts,"TSSundialsSetMinTimeStep_C",(TS,PetscReal),(ts,mindt));CHKERRQ(ierr); 861 PetscFunctionReturn(0); 862 } 863 864 #undef __FUNCT__ 865 #define __FUNCT__ "TSSundialsSetMaxTimeStep" 866 /*@ 867 TSSundialsSetMaxTimeStep - Largest time step to be chosen by the adaptive controller. 868 869 Input Parameter: 870 + ts - the time-step context 871 - maxdt - lowest time step if positive, negative to deactivate 872 873 Level: beginner 874 875 .seealso: TSSundialsSetType(), TSSundialsSetTolerance(), 876 @*/ 877 PetscErrorCode TSSundialsSetMaxTimeStep(TS ts,PetscReal maxdt) 878 { 879 PetscErrorCode ierr; 880 881 PetscFunctionBegin; 882 ierr = PetscTryMethod(ts,"TSSundialsSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));CHKERRQ(ierr); 883 PetscFunctionReturn(0); 884 } 885 886 #undef __FUNCT__ 887 #define __FUNCT__ "TSSundialsMonitorInternalSteps" 888 /*@ 889 TSSundialsMonitorInternalSteps - Monitor Sundials internal steps (Defaults to false). 890 891 Input Parameter: 892 + ts - the time-step context 893 - ft - PETSC_TRUE if monitor, else PETSC_FALSE 894 895 Level: beginner 896 897 .seealso:TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(), 898 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 899 TSSundialsGetIterations(), TSSundialsSetType(), 900 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC() 901 @*/ 902 PetscErrorCode TSSundialsMonitorInternalSteps(TS ts,PetscBool ft) 903 { 904 PetscErrorCode ierr; 905 906 PetscFunctionBegin; 907 ierr = PetscTryMethod(ts,"TSSundialsMonitorInternalSteps_C",(TS,PetscBool),(ts,ft));CHKERRQ(ierr); 908 PetscFunctionReturn(0); 909 } 910 /* -------------------------------------------------------------------------------------------*/ 911 /*MC 912 TSSUNDIALS - ODE solver using the LLNL CVODE/SUNDIALS package (now called SUNDIALS) 913 914 Options Database: 915 + -ts_sundials_type <bdf,adams> 916 . -ts_sundials_gramschmidt_type <modified, classical> - type of orthogonalization inside GMRES 917 . -ts_sundials_atol <tol> - Absolute tolerance for convergence 918 . -ts_sundials_rtol <tol> - Relative tolerance for convergence 919 . -ts_sundials_linear_tolerance <tol> 920 . -ts_sundials_maxl <maxl> - Max dimension of the Krylov subspace 921 - -ts_sundials_monitor_steps - Monitor SUNDIALS internel steps 922 923 924 Notes: This uses its own nonlinear solver and Krylov method so PETSc SNES and KSP options do not apply 925 only PETSc PC options 926 927 Level: beginner 928 929 .seealso: TSCreate(), TS, TSSetType(), TSSundialsSetType(), TSSundialsSetMaxl(), TSSundialsSetLinearTolerance(), 930 TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), TSSundialsGetPC(), TSSundialsGetIterations(), TSSetExactFinalTime() 931 932 M*/ 933 #undef __FUNCT__ 934 #define __FUNCT__ "TSCreate_Sundials" 935 PETSC_EXTERN PetscErrorCode TSCreate_Sundials(TS ts) 936 { 937 TS_Sundials *cvode; 938 PetscErrorCode ierr; 939 PC pc; 940 941 PetscFunctionBegin; 942 ts->ops->reset = TSReset_Sundials; 943 ts->ops->destroy = TSDestroy_Sundials; 944 ts->ops->view = TSView_Sundials; 945 ts->ops->setup = TSSetUp_Sundials; 946 ts->ops->step = TSStep_Sundials; 947 ts->ops->interpolate = TSInterpolate_Sundials; 948 ts->ops->setfromoptions = TSSetFromOptions_Sundials; 949 950 ierr = PetscNewLog(ts,&cvode);CHKERRQ(ierr); 951 952 ts->data = (void*)cvode; 953 cvode->cvode_type = SUNDIALS_BDF; 954 cvode->gtype = SUNDIALS_CLASSICAL_GS; 955 cvode->maxl = 5; 956 cvode->linear_tol = .05; 957 cvode->monitorstep = PETSC_TRUE; 958 959 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)ts),&(cvode->comm_sundials));CHKERRQ(ierr); 960 961 cvode->mindt = -1.; 962 cvode->maxdt = -1.; 963 964 /* set tolerance for Sundials */ 965 cvode->reltol = 1e-6; 966 cvode->abstol = 1e-6; 967 968 /* set PCNONE as default pctype */ 969 ierr = TSSundialsGetPC_Sundials(ts,&pc);CHKERRQ(ierr); 970 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 971 972 if (ts->exact_final_time != TS_EXACTFINALTIME_STEPOVER) ts->exact_final_time = TS_EXACTFINALTIME_STEPOVER; 973 974 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",TSSundialsSetType_Sundials);CHKERRQ(ierr); 975 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",TSSundialsSetMaxl_Sundials);CHKERRQ(ierr); 976 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",TSSundialsSetLinearTolerance_Sundials);CHKERRQ(ierr); 977 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",TSSundialsSetGramSchmidtType_Sundials);CHKERRQ(ierr); 978 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",TSSundialsSetTolerance_Sundials);CHKERRQ(ierr); 979 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",TSSundialsSetMinTimeStep_Sundials);CHKERRQ(ierr); 980 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",TSSundialsSetMaxTimeStep_Sundials);CHKERRQ(ierr); 981 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",TSSundialsGetPC_Sundials);CHKERRQ(ierr); 982 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",TSSundialsGetIterations_Sundials);CHKERRQ(ierr); 983 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",TSSundialsMonitorInternalSteps_Sundials);CHKERRQ(ierr); 984 PetscFunctionReturn(0); 985 } 986