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 nits,lits,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 /* We would like to TSPreStage() and TSPostStage() 137 * before each stage solve but CVode does not appear to support this. */ 138 if (cvode->monitorstep) 139 flag = CVode(mem,tout,cvode->y,&t,CV_ONE_STEP); 140 else 141 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 maxstep %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 /* log inner nonlinear and linear iterations */ 196 ierr = CVodeGetNumNonlinSolvIters(mem,&nits);CHKERRQ(ierr); 197 ierr = CVSpilsGetNumLinIters(mem,&lits);CHKERRQ(ierr); 198 ts->snes_its += nits; ts->ksp_its = lits; 199 200 /* copy the solution from cvode->y to cvode->update and sol */ 201 ierr = VecPlaceArray(cvode->w1,y_data);CHKERRQ(ierr); 202 ierr = VecCopy(cvode->w1,cvode->update);CHKERRQ(ierr); 203 ierr = VecResetArray(cvode->w1);CHKERRQ(ierr); 204 ierr = VecCopy(cvode->update,ts->vec_sol);CHKERRQ(ierr); 205 206 ts->time_step = t - ts->ptime; 207 ts->ptime = t; 208 209 ierr = CVodeGetNumSteps(mem,&nsteps);CHKERRQ(ierr); 210 if (!cvode->monitorstep) ts->steps = nsteps - 1; /* TSStep() increments the step counter by one */ 211 PetscFunctionReturn(0); 212 } 213 214 #undef __FUNCT__ 215 #define __FUNCT__ "TSInterpolate_Sundials" 216 static PetscErrorCode TSInterpolate_Sundials(TS ts,PetscReal t,Vec X) 217 { 218 TS_Sundials *cvode = (TS_Sundials*)ts->data; 219 N_Vector y; 220 PetscErrorCode ierr; 221 PetscScalar *x_data; 222 PetscInt glosize,locsize; 223 224 PetscFunctionBegin; 225 /* get the vector size */ 226 ierr = VecGetSize(X,&glosize);CHKERRQ(ierr); 227 ierr = VecGetLocalSize(X,&locsize);CHKERRQ(ierr); 228 229 /* allocate the memory for N_Vec y */ 230 y = N_VNew_Parallel(cvode->comm_sundials,locsize,glosize); 231 if (!y) SETERRQ(PETSC_COMM_SELF,1,"Interpolated y is not allocated"); 232 233 ierr = VecGetArray(X,&x_data);CHKERRQ(ierr); 234 N_VSetArrayPointer((realtype*)x_data,y); 235 ierr = CVodeGetDky(cvode->mem,t,0,y);CHKERRQ(ierr); 236 ierr = VecRestoreArray(X,&x_data);CHKERRQ(ierr); 237 PetscFunctionReturn(0); 238 } 239 240 #undef __FUNCT__ 241 #define __FUNCT__ "TSReset_Sundials" 242 PetscErrorCode TSReset_Sundials(TS ts) 243 { 244 TS_Sundials *cvode = (TS_Sundials*)ts->data; 245 PetscErrorCode ierr; 246 247 PetscFunctionBegin; 248 ierr = VecDestroy(&cvode->update);CHKERRQ(ierr); 249 ierr = VecDestroy(&cvode->ydot);CHKERRQ(ierr); 250 ierr = VecDestroy(&cvode->w1);CHKERRQ(ierr); 251 ierr = VecDestroy(&cvode->w2);CHKERRQ(ierr); 252 if (cvode->mem) CVodeFree(&cvode->mem); 253 PetscFunctionReturn(0); 254 } 255 256 #undef __FUNCT__ 257 #define __FUNCT__ "TSDestroy_Sundials" 258 PetscErrorCode TSDestroy_Sundials(TS ts) 259 { 260 TS_Sundials *cvode = (TS_Sundials*)ts->data; 261 PetscErrorCode ierr; 262 263 PetscFunctionBegin; 264 ierr = TSReset_Sundials(ts);CHKERRQ(ierr); 265 ierr = MPI_Comm_free(&(cvode->comm_sundials));CHKERRQ(ierr); 266 ierr = PetscFree(ts->data);CHKERRQ(ierr); 267 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",NULL);CHKERRQ(ierr); 268 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",NULL);CHKERRQ(ierr); 269 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",NULL);CHKERRQ(ierr); 270 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",NULL);CHKERRQ(ierr); 271 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",NULL);CHKERRQ(ierr); 272 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",NULL);CHKERRQ(ierr); 273 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",NULL);CHKERRQ(ierr); 274 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",NULL);CHKERRQ(ierr); 275 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",NULL);CHKERRQ(ierr); 276 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",NULL);CHKERRQ(ierr); 277 PetscFunctionReturn(0); 278 } 279 280 #undef __FUNCT__ 281 #define __FUNCT__ "TSSetUp_Sundials" 282 PetscErrorCode TSSetUp_Sundials(TS ts) 283 { 284 TS_Sundials *cvode = (TS_Sundials*)ts->data; 285 PetscErrorCode ierr; 286 PetscInt glosize,locsize,i,flag; 287 PetscScalar *y_data,*parray; 288 void *mem; 289 PC pc; 290 PCType pctype; 291 PetscBool pcnone; 292 293 PetscFunctionBegin; 294 if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for exact final time option 'MATCHSTEP' when using Sundials"); 295 296 /* get the vector size */ 297 ierr = VecGetSize(ts->vec_sol,&glosize);CHKERRQ(ierr); 298 ierr = VecGetLocalSize(ts->vec_sol,&locsize);CHKERRQ(ierr); 299 300 /* allocate the memory for N_Vec y */ 301 cvode->y = N_VNew_Parallel(cvode->comm_sundials,locsize,glosize); 302 if (!cvode->y) SETERRQ(PETSC_COMM_SELF,1,"cvode->y is not allocated"); 303 304 /* initialize N_Vec y: copy ts->vec_sol to cvode->y */ 305 ierr = VecGetArray(ts->vec_sol,&parray);CHKERRQ(ierr); 306 y_data = (PetscScalar*) N_VGetArrayPointer(cvode->y); 307 for (i = 0; i < locsize; i++) y_data[i] = parray[i]; 308 ierr = VecRestoreArray(ts->vec_sol,NULL);CHKERRQ(ierr); 309 310 ierr = VecDuplicate(ts->vec_sol,&cvode->update);CHKERRQ(ierr); 311 ierr = VecDuplicate(ts->vec_sol,&cvode->ydot);CHKERRQ(ierr); 312 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->update);CHKERRQ(ierr); 313 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->ydot);CHKERRQ(ierr); 314 315 /* 316 Create work vectors for the TSPSolve_Sundials() routine. Note these are 317 allocated with zero space arrays because the actual array space is provided 318 by Sundials and set using VecPlaceArray(). 319 */ 320 ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,0,&cvode->w1);CHKERRQ(ierr); 321 ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,0,&cvode->w2);CHKERRQ(ierr); 322 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w1);CHKERRQ(ierr); 323 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w2);CHKERRQ(ierr); 324 325 /* Call CVodeCreate to create the solver memory and the use of a Newton iteration */ 326 mem = CVodeCreate(cvode->cvode_type, CV_NEWTON); 327 if (!mem) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"CVodeCreate() fails"); 328 cvode->mem = mem; 329 330 /* Set the pointer to user-defined data */ 331 flag = CVodeSetUserData(mem, ts); 332 if (flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSetUserData() fails"); 333 334 /* Sundials may choose to use a smaller initial step, but will never use a larger step. */ 335 flag = CVodeSetInitStep(mem,(realtype)ts->time_step); 336 if (flag) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetInitStep() failed"); 337 if (cvode->mindt > 0) { 338 flag = CVodeSetMinStep(mem,(realtype)cvode->mindt); 339 if (flag) { 340 if (flag == CV_MEM_NULL) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed, cvode_mem pointer is NULL"); 341 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"); 342 else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed"); 343 } 344 } 345 if (cvode->maxdt > 0) { 346 flag = CVodeSetMaxStep(mem,(realtype)cvode->maxdt); 347 if (flag) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMaxStep() failed"); 348 } 349 350 /* Call CVodeInit to initialize the integrator memory and specify the 351 * user's right hand side function in u'=f(t,u), the inital time T0, and 352 * the initial dependent variable vector cvode->y */ 353 flag = CVodeInit(mem,TSFunction_Sundials,ts->ptime,cvode->y); 354 if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeInit() fails, flag %d",flag); 355 356 /* specifies scalar relative and absolute tolerances */ 357 flag = CVodeSStolerances(mem,cvode->reltol,cvode->abstol); 358 if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSStolerances() fails, flag %d",flag); 359 360 /* Specify max num of steps to be taken by cvode in its attempt to reach the next output time */ 361 flag = CVodeSetMaxNumSteps(mem,ts->max_steps); 362 if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSetMaxNumSteps() fails, flag %d",flag); 363 364 /* call CVSpgmr to use GMRES as the linear solver. */ 365 /* setup the ode integrator with the given preconditioner */ 366 ierr = TSSundialsGetPC(ts,&pc);CHKERRQ(ierr); 367 ierr = PCGetType(pc,&pctype);CHKERRQ(ierr); 368 ierr = PetscObjectTypeCompare((PetscObject)pc,PCNONE,&pcnone);CHKERRQ(ierr); 369 if (pcnone) { 370 flag = CVSpgmr(mem,PREC_NONE,0); 371 if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmr() fails, flag %d",flag); 372 } else { 373 flag = CVSpgmr(mem,PREC_LEFT,cvode->maxl); 374 if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmr() fails, flag %d",flag); 375 376 /* Set preconditioner and solve routines Precond and PSolve, 377 and the pointer to the user-defined block data */ 378 flag = CVSpilsSetPreconditioner(mem,TSPrecond_Sundials,TSPSolve_Sundials); 379 if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpilsSetPreconditioner() fails, flag %d", flag); 380 } 381 382 flag = CVSpilsSetGSType(mem, MODIFIED_GS); 383 if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmrSetGSType() fails, flag %d",flag); 384 PetscFunctionReturn(0); 385 } 386 387 /* type of CVODE linear multistep method */ 388 const char *const TSSundialsLmmTypes[] = {"","ADAMS","BDF","TSSundialsLmmType","SUNDIALS_",0}; 389 /* type of G-S orthogonalization used by CVODE linear solver */ 390 const char *const TSSundialsGramSchmidtTypes[] = {"","MODIFIED","CLASSICAL","TSSundialsGramSchmidtType","SUNDIALS_",0}; 391 392 #undef __FUNCT__ 393 #define __FUNCT__ "TSSetFromOptions_Sundials" 394 PetscErrorCode TSSetFromOptions_Sundials(PetscOptionItems *PetscOptionsObject,TS ts) 395 { 396 TS_Sundials *cvode = (TS_Sundials*)ts->data; 397 PetscErrorCode ierr; 398 int indx; 399 PetscBool flag; 400 PC pc; 401 402 PetscFunctionBegin; 403 ierr = PetscOptionsHead(PetscOptionsObject,"SUNDIALS ODE solver options");CHKERRQ(ierr); 404 ierr = PetscOptionsEList("-ts_sundials_type","Scheme","TSSundialsSetType",TSSundialsLmmTypes,3,TSSundialsLmmTypes[cvode->cvode_type],&indx,&flag);CHKERRQ(ierr); 405 if (flag) { 406 ierr = TSSundialsSetType(ts,(TSSundialsLmmType)indx);CHKERRQ(ierr); 407 } 408 ierr = PetscOptionsEList("-ts_sundials_gramschmidt_type","Type of orthogonalization","TSSundialsSetGramSchmidtType",TSSundialsGramSchmidtTypes,3,TSSundialsGramSchmidtTypes[cvode->gtype],&indx,&flag);CHKERRQ(ierr); 409 if (flag) { 410 ierr = TSSundialsSetGramSchmidtType(ts,(TSSundialsGramSchmidtType)indx);CHKERRQ(ierr); 411 } 412 ierr = PetscOptionsReal("-ts_sundials_atol","Absolute tolerance for convergence","TSSundialsSetTolerance",cvode->abstol,&cvode->abstol,NULL);CHKERRQ(ierr); 413 ierr = PetscOptionsReal("-ts_sundials_rtol","Relative tolerance for convergence","TSSundialsSetTolerance",cvode->reltol,&cvode->reltol,NULL);CHKERRQ(ierr); 414 ierr = PetscOptionsReal("-ts_sundials_mindt","Minimum step size","TSSundialsSetMinTimeStep",cvode->mindt,&cvode->mindt,NULL);CHKERRQ(ierr); 415 ierr = PetscOptionsReal("-ts_sundials_maxdt","Maximum step size","TSSundialsSetMaxTimeStep",cvode->maxdt,&cvode->maxdt,NULL);CHKERRQ(ierr); 416 ierr = PetscOptionsReal("-ts_sundials_linear_tolerance","Convergence tolerance for linear solve","TSSundialsSetLinearTolerance",cvode->linear_tol,&cvode->linear_tol,NULL);CHKERRQ(ierr); 417 ierr = PetscOptionsInt("-ts_sundials_maxl","Max dimension of the Krylov subspace","TSSundialsSetMaxl",cvode->maxl,&cvode->maxl,NULL);CHKERRQ(ierr); 418 ierr = PetscOptionsBool("-ts_sundials_monitor_steps","Monitor SUNDIALS internal steps","TSSundialsMonitorInternalSteps",cvode->monitorstep,&cvode->monitorstep,NULL);CHKERRQ(ierr); 419 ierr = PetscOptionsTail();CHKERRQ(ierr); 420 ierr = TSSundialsGetPC(ts,&pc);CHKERRQ(ierr); 421 ierr = PCSetFromOptions(pc);CHKERRQ(ierr); 422 PetscFunctionReturn(0); 423 } 424 425 #undef __FUNCT__ 426 #define __FUNCT__ "TSView_Sundials" 427 PetscErrorCode TSView_Sundials(TS ts,PetscViewer viewer) 428 { 429 TS_Sundials *cvode = (TS_Sundials*)ts->data; 430 PetscErrorCode ierr; 431 char *type; 432 char atype[] = "Adams"; 433 char btype[] = "BDF: backward differentiation formula"; 434 PetscBool iascii,isstring; 435 long int nsteps,its,nfevals,nlinsetups,nfails,itmp; 436 PetscInt qlast,qcur; 437 PetscReal hinused,hlast,hcur,tcur,tolsfac; 438 PC pc; 439 440 PetscFunctionBegin; 441 if (cvode->cvode_type == SUNDIALS_ADAMS) type = atype; 442 else type = btype; 443 444 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 445 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 446 if (iascii) { 447 ierr = PetscViewerASCIIPrintf(viewer,"Sundials integrater does not use SNES!\n");CHKERRQ(ierr); 448 ierr = PetscViewerASCIIPrintf(viewer,"Sundials integrater type %s\n",type);CHKERRQ(ierr); 449 ierr = PetscViewerASCIIPrintf(viewer,"Sundials abs tol %g rel tol %g\n",cvode->abstol,cvode->reltol);CHKERRQ(ierr); 450 ierr = PetscViewerASCIIPrintf(viewer,"Sundials linear solver tolerance factor %g\n",cvode->linear_tol);CHKERRQ(ierr); 451 ierr = PetscViewerASCIIPrintf(viewer,"Sundials max dimension of Krylov subspace %D\n",cvode->maxl);CHKERRQ(ierr); 452 if (cvode->gtype == SUNDIALS_MODIFIED_GS) { 453 ierr = PetscViewerASCIIPrintf(viewer,"Sundials using modified Gram-Schmidt for orthogonalization in GMRES\n");CHKERRQ(ierr); 454 } else { 455 ierr = PetscViewerASCIIPrintf(viewer,"Sundials using unmodified (classical) Gram-Schmidt for orthogonalization in GMRES\n");CHKERRQ(ierr); 456 } 457 if (cvode->mindt > 0) {ierr = PetscViewerASCIIPrintf(viewer,"Sundials minimum time step %g\n",cvode->mindt);CHKERRQ(ierr);} 458 if (cvode->maxdt > 0) {ierr = PetscViewerASCIIPrintf(viewer,"Sundials maximum time step %g\n",cvode->maxdt);CHKERRQ(ierr);} 459 460 /* Outputs from CVODE, CVSPILS */ 461 ierr = CVodeGetTolScaleFactor(cvode->mem,&tolsfac);CHKERRQ(ierr); 462 ierr = PetscViewerASCIIPrintf(viewer,"Sundials suggested factor for tolerance scaling %g\n",tolsfac);CHKERRQ(ierr); 463 ierr = CVodeGetIntegratorStats(cvode->mem,&nsteps,&nfevals, 464 &nlinsetups,&nfails,&qlast,&qcur, 465 &hinused,&hlast,&hcur,&tcur);CHKERRQ(ierr); 466 ierr = PetscViewerASCIIPrintf(viewer,"Sundials cumulative number of internal steps %D\n",nsteps);CHKERRQ(ierr); 467 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to rhs function %D\n",nfevals);CHKERRQ(ierr); 468 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to linear solver setup function %D\n",nlinsetups);CHKERRQ(ierr); 469 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of error test failures %D\n",nfails);CHKERRQ(ierr); 470 471 ierr = CVodeGetNonlinSolvStats(cvode->mem,&its,&nfails);CHKERRQ(ierr); 472 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear solver iterations %D\n",its);CHKERRQ(ierr); 473 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear convergence failure %D\n",nfails);CHKERRQ(ierr); 474 475 ierr = CVSpilsGetNumLinIters(cvode->mem, &its);CHKERRQ(ierr); /* its = no. of calls to TSPrecond_Sundials() */ 476 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of linear iterations %D\n",its);CHKERRQ(ierr); 477 ierr = CVSpilsGetNumConvFails(cvode->mem,&itmp);CHKERRQ(ierr); 478 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of linear convergence failures %D\n",itmp);CHKERRQ(ierr); 479 480 ierr = TSSundialsGetPC(ts,&pc);CHKERRQ(ierr); 481 ierr = PCView(pc,viewer);CHKERRQ(ierr); 482 ierr = CVSpilsGetNumPrecEvals(cvode->mem,&itmp);CHKERRQ(ierr); 483 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner evaluations %D\n",itmp);CHKERRQ(ierr); 484 ierr = CVSpilsGetNumPrecSolves(cvode->mem,&itmp);CHKERRQ(ierr); 485 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner solves %D\n",itmp);CHKERRQ(ierr); 486 487 ierr = CVSpilsGetNumJtimesEvals(cvode->mem,&itmp);CHKERRQ(ierr); 488 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of Jacobian-vector product evaluations %D\n",itmp);CHKERRQ(ierr); 489 ierr = CVSpilsGetNumRhsEvals(cvode->mem,&itmp);CHKERRQ(ierr); 490 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of rhs calls for finite diff. Jacobian-vector evals %D\n",itmp);CHKERRQ(ierr); 491 } else if (isstring) { 492 ierr = PetscViewerStringSPrintf(viewer,"Sundials type %s",type);CHKERRQ(ierr); 493 } 494 PetscFunctionReturn(0); 495 } 496 497 498 /* --------------------------------------------------------------------------*/ 499 #undef __FUNCT__ 500 #define __FUNCT__ "TSSundialsSetType_Sundials" 501 PetscErrorCode TSSundialsSetType_Sundials(TS ts,TSSundialsLmmType type) 502 { 503 TS_Sundials *cvode = (TS_Sundials*)ts->data; 504 505 PetscFunctionBegin; 506 cvode->cvode_type = type; 507 PetscFunctionReturn(0); 508 } 509 510 #undef __FUNCT__ 511 #define __FUNCT__ "TSSundialsSetMaxl_Sundials" 512 PetscErrorCode TSSundialsSetMaxl_Sundials(TS ts,PetscInt maxl) 513 { 514 TS_Sundials *cvode = (TS_Sundials*)ts->data; 515 516 PetscFunctionBegin; 517 cvode->maxl = maxl; 518 PetscFunctionReturn(0); 519 } 520 521 #undef __FUNCT__ 522 #define __FUNCT__ "TSSundialsSetLinearTolerance_Sundials" 523 PetscErrorCode TSSundialsSetLinearTolerance_Sundials(TS ts,double tol) 524 { 525 TS_Sundials *cvode = (TS_Sundials*)ts->data; 526 527 PetscFunctionBegin; 528 cvode->linear_tol = tol; 529 PetscFunctionReturn(0); 530 } 531 532 #undef __FUNCT__ 533 #define __FUNCT__ "TSSundialsSetGramSchmidtType_Sundials" 534 PetscErrorCode TSSundialsSetGramSchmidtType_Sundials(TS ts,TSSundialsGramSchmidtType type) 535 { 536 TS_Sundials *cvode = (TS_Sundials*)ts->data; 537 538 PetscFunctionBegin; 539 cvode->gtype = type; 540 PetscFunctionReturn(0); 541 } 542 543 #undef __FUNCT__ 544 #define __FUNCT__ "TSSundialsSetTolerance_Sundials" 545 PetscErrorCode TSSundialsSetTolerance_Sundials(TS ts,double aabs,double rel) 546 { 547 TS_Sundials *cvode = (TS_Sundials*)ts->data; 548 549 PetscFunctionBegin; 550 if (aabs != PETSC_DECIDE) cvode->abstol = aabs; 551 if (rel != PETSC_DECIDE) cvode->reltol = rel; 552 PetscFunctionReturn(0); 553 } 554 555 #undef __FUNCT__ 556 #define __FUNCT__ "TSSundialsSetMinTimeStep_Sundials" 557 PetscErrorCode TSSundialsSetMinTimeStep_Sundials(TS ts,PetscReal mindt) 558 { 559 TS_Sundials *cvode = (TS_Sundials*)ts->data; 560 561 PetscFunctionBegin; 562 cvode->mindt = mindt; 563 PetscFunctionReturn(0); 564 } 565 566 #undef __FUNCT__ 567 #define __FUNCT__ "TSSundialsSetMaxTimeStep_Sundials" 568 PetscErrorCode TSSundialsSetMaxTimeStep_Sundials(TS ts,PetscReal maxdt) 569 { 570 TS_Sundials *cvode = (TS_Sundials*)ts->data; 571 572 PetscFunctionBegin; 573 cvode->maxdt = maxdt; 574 PetscFunctionReturn(0); 575 } 576 #undef __FUNCT__ 577 #define __FUNCT__ "TSSundialsGetPC_Sundials" 578 PetscErrorCode TSSundialsGetPC_Sundials(TS ts,PC *pc) 579 { 580 SNES snes; 581 KSP ksp; 582 PetscErrorCode ierr; 583 584 PetscFunctionBegin; 585 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 586 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 587 ierr = KSPGetPC(ksp,pc);CHKERRQ(ierr); 588 PetscFunctionReturn(0); 589 } 590 591 #undef __FUNCT__ 592 #define __FUNCT__ "TSSundialsGetIterations_Sundials" 593 PetscErrorCode TSSundialsGetIterations_Sundials(TS ts,int *nonlin,int *lin) 594 { 595 PetscFunctionBegin; 596 if (nonlin) *nonlin = ts->snes_its; 597 if (lin) *lin = ts->ksp_its; 598 PetscFunctionReturn(0); 599 } 600 601 #undef __FUNCT__ 602 #define __FUNCT__ "TSSundialsMonitorInternalSteps_Sundials" 603 PetscErrorCode TSSundialsMonitorInternalSteps_Sundials(TS ts,PetscBool s) 604 { 605 TS_Sundials *cvode = (TS_Sundials*)ts->data; 606 607 PetscFunctionBegin; 608 cvode->monitorstep = s; 609 PetscFunctionReturn(0); 610 } 611 /* -------------------------------------------------------------------------------------------*/ 612 613 #undef __FUNCT__ 614 #define __FUNCT__ "TSSundialsGetIterations" 615 /*@C 616 TSSundialsGetIterations - Gets the number of nonlinear and linear iterations used so far by Sundials. 617 618 Not Collective 619 620 Input parameters: 621 . ts - the time-step context 622 623 Output Parameters: 624 + nonlin - number of nonlinear iterations 625 - lin - number of linear iterations 626 627 Level: advanced 628 629 Notes: 630 These return the number since the creation of the TS object 631 632 .keywords: non-linear iterations, linear iterations 633 634 .seealso: TSSundialsSetType(), TSSundialsSetMaxl(), 635 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 636 TSSundialsGetIterations(), TSSundialsSetType(), 637 TSSundialsSetLinearTolerance(), TSSundialsGetPC(), TSSetExactFinalTime() 638 639 @*/ 640 PetscErrorCode TSSundialsGetIterations(TS ts,int *nonlin,int *lin) 641 { 642 PetscErrorCode ierr; 643 644 PetscFunctionBegin; 645 ierr = PetscUseMethod(ts,"TSSundialsGetIterations_C",(TS,int*,int*),(ts,nonlin,lin));CHKERRQ(ierr); 646 PetscFunctionReturn(0); 647 } 648 649 #undef __FUNCT__ 650 #define __FUNCT__ "TSSundialsSetType" 651 /*@ 652 TSSundialsSetType - Sets the method that Sundials will use for integration. 653 654 Logically Collective on TS 655 656 Input parameters: 657 + ts - the time-step context 658 - type - one of SUNDIALS_ADAMS or SUNDIALS_BDF 659 660 Level: intermediate 661 662 .keywords: Adams, backward differentiation formula 663 664 .seealso: TSSundialsGetIterations(), TSSundialsSetMaxl(), 665 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 666 TSSundialsGetIterations(), TSSundialsSetType(), 667 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(), 668 TSSetExactFinalTime() 669 @*/ 670 PetscErrorCode TSSundialsSetType(TS ts,TSSundialsLmmType type) 671 { 672 PetscErrorCode ierr; 673 674 PetscFunctionBegin; 675 ierr = PetscTryMethod(ts,"TSSundialsSetType_C",(TS,TSSundialsLmmType),(ts,type));CHKERRQ(ierr); 676 PetscFunctionReturn(0); 677 } 678 679 #undef __FUNCT__ 680 #define __FUNCT__ "TSSundialsSetMaxl" 681 /*@ 682 TSSundialsSetMaxl - Sets the dimension of the Krylov space used by 683 GMRES in the linear solver in SUNDIALS. SUNDIALS DOES NOT use restarted GMRES so 684 this is the maximum number of GMRES steps that will be used. 685 686 Logically Collective on TS 687 688 Input parameters: 689 + ts - the time-step context 690 - maxl - number of direction vectors (the dimension of Krylov subspace). 691 692 Level: advanced 693 694 .keywords: GMRES 695 696 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), 697 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 698 TSSundialsGetIterations(), TSSundialsSetType(), 699 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(), 700 TSSetExactFinalTime() 701 702 @*/ 703 PetscErrorCode TSSundialsSetMaxl(TS ts,PetscInt maxl) 704 { 705 PetscErrorCode ierr; 706 707 PetscFunctionBegin; 708 PetscValidLogicalCollectiveInt(ts,maxl,2); 709 ierr = PetscTryMethod(ts,"TSSundialsSetMaxl_C",(TS,PetscInt),(ts,maxl));CHKERRQ(ierr); 710 PetscFunctionReturn(0); 711 } 712 713 #undef __FUNCT__ 714 #define __FUNCT__ "TSSundialsSetLinearTolerance" 715 /*@ 716 TSSundialsSetLinearTolerance - Sets the tolerance used to solve the linear 717 system by SUNDIALS. 718 719 Logically Collective on TS 720 721 Input parameters: 722 + ts - the time-step context 723 - tol - the factor by which the tolerance on the nonlinear solver is 724 multiplied to get the tolerance on the linear solver, .05 by default. 725 726 Level: advanced 727 728 .keywords: GMRES, linear convergence tolerance, SUNDIALS 729 730 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(), 731 TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 732 TSSundialsGetIterations(), TSSundialsSetType(), 733 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(), 734 TSSetExactFinalTime() 735 736 @*/ 737 PetscErrorCode TSSundialsSetLinearTolerance(TS ts,double tol) 738 { 739 PetscErrorCode ierr; 740 741 PetscFunctionBegin; 742 PetscValidLogicalCollectiveReal(ts,tol,2); 743 ierr = PetscTryMethod(ts,"TSSundialsSetLinearTolerance_C",(TS,double),(ts,tol));CHKERRQ(ierr); 744 PetscFunctionReturn(0); 745 } 746 747 #undef __FUNCT__ 748 #define __FUNCT__ "TSSundialsSetGramSchmidtType" 749 /*@ 750 TSSundialsSetGramSchmidtType - Sets type of orthogonalization used 751 in GMRES method by SUNDIALS linear solver. 752 753 Logically Collective on TS 754 755 Input parameters: 756 + ts - the time-step context 757 - type - either SUNDIALS_MODIFIED_GS or SUNDIALS_CLASSICAL_GS 758 759 Level: advanced 760 761 .keywords: Sundials, orthogonalization 762 763 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(), 764 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), 765 TSSundialsGetIterations(), TSSundialsSetType(), 766 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(), 767 TSSetExactFinalTime() 768 769 @*/ 770 PetscErrorCode TSSundialsSetGramSchmidtType(TS ts,TSSundialsGramSchmidtType type) 771 { 772 PetscErrorCode ierr; 773 774 PetscFunctionBegin; 775 ierr = PetscTryMethod(ts,"TSSundialsSetGramSchmidtType_C",(TS,TSSundialsGramSchmidtType),(ts,type));CHKERRQ(ierr); 776 PetscFunctionReturn(0); 777 } 778 779 #undef __FUNCT__ 780 #define __FUNCT__ "TSSundialsSetTolerance" 781 /*@ 782 TSSundialsSetTolerance - Sets the absolute and relative tolerance used by 783 Sundials for error control. 784 785 Logically Collective on TS 786 787 Input parameters: 788 + ts - the time-step context 789 . aabs - the absolute tolerance 790 - rel - the relative tolerance 791 792 See the Cvode/Sundials users manual for exact details on these parameters. Essentially 793 these regulate the size of the error for a SINGLE timestep. 794 795 Level: intermediate 796 797 .keywords: Sundials, tolerance 798 799 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESMaxl(), 800 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), 801 TSSundialsGetIterations(), TSSundialsSetType(), 802 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(), 803 TSSetExactFinalTime() 804 805 @*/ 806 PetscErrorCode TSSundialsSetTolerance(TS ts,double aabs,double rel) 807 { 808 PetscErrorCode ierr; 809 810 PetscFunctionBegin; 811 ierr = PetscTryMethod(ts,"TSSundialsSetTolerance_C",(TS,double,double),(ts,aabs,rel));CHKERRQ(ierr); 812 PetscFunctionReturn(0); 813 } 814 815 #undef __FUNCT__ 816 #define __FUNCT__ "TSSundialsGetPC" 817 /*@ 818 TSSundialsGetPC - Extract the PC context from a time-step context for Sundials. 819 820 Input Parameter: 821 . ts - the time-step context 822 823 Output Parameter: 824 . pc - the preconditioner context 825 826 Level: advanced 827 828 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(), 829 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 830 TSSundialsGetIterations(), TSSundialsSetType(), 831 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance() 832 @*/ 833 PetscErrorCode TSSundialsGetPC(TS ts,PC *pc) 834 { 835 PetscErrorCode ierr; 836 837 PetscFunctionBegin; 838 ierr = PetscUseMethod(ts,"TSSundialsGetPC_C",(TS,PC*),(ts,pc));CHKERRQ(ierr); 839 PetscFunctionReturn(0); 840 } 841 842 #undef __FUNCT__ 843 #define __FUNCT__ "TSSundialsSetMinTimeStep" 844 /*@ 845 TSSundialsSetMinTimeStep - Smallest time step to be chosen by the adaptive controller. 846 847 Input Parameter: 848 + ts - the time-step context 849 - mindt - lowest time step if positive, negative to deactivate 850 851 Note: 852 Sundials will error if it is not possible to keep the estimated truncation error below 853 the tolerance set with TSSundialsSetTolerance() without going below this step size. 854 855 Level: beginner 856 857 .seealso: TSSundialsSetType(), TSSundialsSetTolerance(), 858 @*/ 859 PetscErrorCode TSSundialsSetMinTimeStep(TS ts,PetscReal mindt) 860 { 861 PetscErrorCode ierr; 862 863 PetscFunctionBegin; 864 ierr = PetscTryMethod(ts,"TSSundialsSetMinTimeStep_C",(TS,PetscReal),(ts,mindt));CHKERRQ(ierr); 865 PetscFunctionReturn(0); 866 } 867 868 #undef __FUNCT__ 869 #define __FUNCT__ "TSSundialsSetMaxTimeStep" 870 /*@ 871 TSSundialsSetMaxTimeStep - Largest time step to be chosen by the adaptive controller. 872 873 Input Parameter: 874 + ts - the time-step context 875 - maxdt - lowest time step if positive, negative to deactivate 876 877 Level: beginner 878 879 .seealso: TSSundialsSetType(), TSSundialsSetTolerance(), 880 @*/ 881 PetscErrorCode TSSundialsSetMaxTimeStep(TS ts,PetscReal maxdt) 882 { 883 PetscErrorCode ierr; 884 885 PetscFunctionBegin; 886 ierr = PetscTryMethod(ts,"TSSundialsSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));CHKERRQ(ierr); 887 PetscFunctionReturn(0); 888 } 889 890 #undef __FUNCT__ 891 #define __FUNCT__ "TSSundialsMonitorInternalSteps" 892 /*@ 893 TSSundialsMonitorInternalSteps - Monitor Sundials internal steps (Defaults to false). 894 895 Input Parameter: 896 + ts - the time-step context 897 - ft - PETSC_TRUE if monitor, else PETSC_FALSE 898 899 Level: beginner 900 901 .seealso:TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(), 902 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), 903 TSSundialsGetIterations(), TSSundialsSetType(), 904 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC() 905 @*/ 906 PetscErrorCode TSSundialsMonitorInternalSteps(TS ts,PetscBool ft) 907 { 908 PetscErrorCode ierr; 909 910 PetscFunctionBegin; 911 ierr = PetscTryMethod(ts,"TSSundialsMonitorInternalSteps_C",(TS,PetscBool),(ts,ft));CHKERRQ(ierr); 912 PetscFunctionReturn(0); 913 } 914 /* -------------------------------------------------------------------------------------------*/ 915 /*MC 916 TSSUNDIALS - ODE solver using the LLNL CVODE/SUNDIALS package (now called SUNDIALS) 917 918 Options Database: 919 + -ts_sundials_type <bdf,adams> - 920 . -ts_sundials_gramschmidt_type <modified, classical> - type of orthogonalization inside GMRES 921 . -ts_sundials_atol <tol> - Absolute tolerance for convergence 922 . -ts_sundials_rtol <tol> - Relative tolerance for convergence 923 . -ts_sundials_linear_tolerance <tol> - 924 . -ts_sundials_maxl <maxl> - Max dimension of the Krylov subspace 925 - -ts_sundials_monitor_steps - Monitor SUNDIALS internal steps 926 927 928 Notes: This uses its own nonlinear solver and Krylov method so PETSc SNES and KSP options do not apply, 929 only PETSc PC options. 930 931 Level: beginner 932 933 .seealso: TSCreate(), TS, TSSetType(), TSSundialsSetType(), TSSundialsSetMaxl(), TSSundialsSetLinearTolerance(), 934 TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), TSSundialsGetPC(), TSSundialsGetIterations(), TSSetExactFinalTime() 935 936 M*/ 937 #undef __FUNCT__ 938 #define __FUNCT__ "TSCreate_Sundials" 939 PETSC_EXTERN PetscErrorCode TSCreate_Sundials(TS ts) 940 { 941 TS_Sundials *cvode; 942 PetscErrorCode ierr; 943 PC pc; 944 945 PetscFunctionBegin; 946 ts->ops->reset = TSReset_Sundials; 947 ts->ops->destroy = TSDestroy_Sundials; 948 ts->ops->view = TSView_Sundials; 949 ts->ops->setup = TSSetUp_Sundials; 950 ts->ops->step = TSStep_Sundials; 951 ts->ops->interpolate = TSInterpolate_Sundials; 952 ts->ops->setfromoptions = TSSetFromOptions_Sundials; 953 954 ierr = PetscNewLog(ts,&cvode);CHKERRQ(ierr); 955 956 ts->data = (void*)cvode; 957 cvode->cvode_type = SUNDIALS_BDF; 958 cvode->gtype = SUNDIALS_CLASSICAL_GS; 959 cvode->maxl = 5; 960 cvode->linear_tol = .05; 961 cvode->monitorstep = PETSC_TRUE; 962 963 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)ts),&(cvode->comm_sundials));CHKERRQ(ierr); 964 965 cvode->mindt = -1.; 966 cvode->maxdt = -1.; 967 968 /* set tolerance for Sundials */ 969 cvode->reltol = 1e-6; 970 cvode->abstol = 1e-6; 971 972 /* set PCNONE as default pctype */ 973 ierr = TSSundialsGetPC_Sundials(ts,&pc);CHKERRQ(ierr); 974 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 975 976 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",TSSundialsSetType_Sundials);CHKERRQ(ierr); 977 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",TSSundialsSetMaxl_Sundials);CHKERRQ(ierr); 978 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",TSSundialsSetLinearTolerance_Sundials);CHKERRQ(ierr); 979 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",TSSundialsSetGramSchmidtType_Sundials);CHKERRQ(ierr); 980 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",TSSundialsSetTolerance_Sundials);CHKERRQ(ierr); 981 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",TSSundialsSetMinTimeStep_Sundials);CHKERRQ(ierr); 982 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",TSSundialsSetMaxTimeStep_Sundials);CHKERRQ(ierr); 983 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",TSSundialsGetPC_Sundials);CHKERRQ(ierr); 984 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",TSSundialsGetIterations_Sundials);CHKERRQ(ierr); 985 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",TSSundialsMonitorInternalSteps_Sundials);CHKERRQ(ierr); 986 PetscFunctionReturn(0); 987 } 988