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