Actual source code: petscpvode.c

      1: #define PETSCTS_DLL

      3: /*
      4:     Provides a PETSc interface to PVODE. Alan Hindmarsh's parallel ODE
      5:    solver.
      6: */

      8: #include "src/ts/impls/implicit/pvode/petscpvode.h"  /*I "petscts.h" I*/

     12: /*@C
     13:    TSPVodeGetParameters - Extract "iopt" and "ropt" PVODE parameters

     15:    Input Parameter:
     16: .    ts - the time-step context

     18:    Output Parameters:
     19: +   opt_size - size of the parameter arrays (will be set to PVODE value OPT_SIZE)
     20: .   iopt - the integer parameters
     21: -   ropt - the double paramters

     24:    Level: advanced
     25: g
     26:    Notes: You may pass PETSC_NULL for any value you do not desire

     28:           PETSc initializes these array with the default PVODE values of 0, you may
     29:           change them before calling the TS solver.

     31:           See the PVODE include file cvode.h for the definitions of the fields

     33:     Suggested by: Timothy William Chevalier

     35: .seealso: TSPVodeGetIterations(), TSPVodeSetType(), TSPVodeSetGMRESRestart(),
     36:           TSPVodeSetLinearTolerance(), TSPVodeSetGramSchmidtType(), TSPVodeSetTolerance(),
     37:           TSPVodeGetIterations(), TSPVodeSetType(), TSPVodeSetGMRESRestart(),
     38:           TSPVodeSetLinearTolerance(), TSPVodeSetTolerance()
     39: @*/
     40: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeGetParameters(TS ts,int *opt_size,long int *iopt[],double *ropt[])
     41: {
     42:   TS_PVode     *cvode = (TS_PVode*)ts->data;

     45:   if (opt_size) *opt_size = OPT_SIZE;
     46:   if (iopt)     *iopt     = cvode->iopt;
     47:   if (ropt)     *ropt     = cvode->ropt;
     48:   return(0);
     49: }

     51: /*
     52:       TSPrecond_PVode - function that we provide to PVODE to
     53:                         evaluate the preconditioner.

     55:     Contributed by: Liyang Xu

     57: */
     60: PetscErrorCode TSPrecond_PVode(integertype N,realtype tn,N_Vector y,N_Vector fy,booleantype jok,
     61:                     booleantype *jcurPtr,realtype _gamma,N_Vector ewt,realtype h,
     62:                     realtype uround,long int *nfePtr,void *P_data,
     63:                     N_Vector vtemp1,N_Vector vtemp2,N_Vector vtemp3)
     64: {
     65:   TS           ts = (TS) P_data;
     66:   TS_PVode     *cvode = (TS_PVode*)ts->data;
     67:   PC           pc = cvode->pc;
     69:   Mat          Jac = ts->B;
     70:   Vec          tmpy = cvode->w1;
     71:   PetscScalar  one = 1.0,gm;
     72:   MatStructure str = DIFFERENT_NONZERO_PATTERN;
     73:
     75:   /* This allows us to construct preconditioners in-place if we like */
     76:   MatSetUnfactored(Jac);

     78:   /*
     79:        jok - TRUE means reuse current Jacobian else recompute Jacobian
     80:   */
     81:   if (jok) {
     82:     MatCopy(cvode->pmat,Jac,SAME_NONZERO_PATTERN);
     83:     str      = SAME_NONZERO_PATTERN;
     84:     *jcurPtr = FALSE;
     85:   } else {
     86:     /* make PETSc vector tmpy point to PVODE vector y */
     87:     VecPlaceArray(tmpy,N_VGetData(y));

     89:     /* compute the Jacobian */
     90:     TSComputeRHSJacobian(ts,ts->ptime,tmpy,&Jac,&Jac,&str);

     92:     /* copy the Jacobian matrix */
     93:     if (!cvode->pmat) {
     94:       MatDuplicate(Jac,MAT_COPY_VALUES,&cvode->pmat);
     95:       PetscLogObjectParent(ts,cvode->pmat);
     96:     }
     97:     MatCopy(Jac,cvode->pmat,SAME_NONZERO_PATTERN);

     99:     *jcurPtr = TRUE;
    100:   }

    102:   /* construct I-gamma*Jac  */
    103:   gm   = -_gamma;
    104:   MatScale(Jac,gm);
    105:   MatShift(Jac,one);
    106:
    107:   PCSetOperators(pc,Jac,Jac,str);
    108:   return(0);
    109: }

    111: /*
    112:      TSPSolve_PVode -  routine that we provide to PVode that applies the preconditioner.
    113:
    114:     Contributed by: Liyang Xu

    116: */
    119: PetscErrorCode TSPSolve_PVode(integertype N,realtype tn,N_Vector y,N_Vector fy,N_Vector vtemp,
    120:                    realtype _gamma,N_Vector ewt,realtype delta,long int *nfePtr,
    121:                    N_Vector r,int lr,void *P_data,N_Vector z)
    122: {
    123:   TS       ts = (TS) P_data;
    124:   TS_PVode *cvode = (TS_PVode*)ts->data;
    125:   PC       pc = cvode->pc;
    126:   Vec      rr = cvode->w1,xx = cvode->w2;

    130:   /*
    131:       Make the PETSc work vectors rr and xx point to the arrays in the PVODE vectors
    132:   */
    133:   VecPlaceArray(rr,N_VGetData(r));
    134:   VecPlaceArray(xx,N_VGetData(z));

    136:   /*
    137:       Solve the Px=r and put the result in xx
    138:   */
    139:   PCApply(pc,rr,xx);
    140:   cvode->linear_solves++;

    143:   return(0);
    144: }

    146: /*
    147:         TSFunction_PVode - routine that we provide to PVode that applies the right hand side.
    148:
    149:     Contributed by: Liyang Xu
    150: */
    153: void TSFunction_PVode(int N,double t,N_Vector y,N_Vector ydot,void *ctx)
    154: {
    155:   TS        ts = (TS) ctx;
    156:   TS_PVode *cvode = (TS_PVode*)ts->data;
    157:   Vec       tmpx = cvode->w1,tmpy = cvode->w2;

    161:   /*
    162:       Make the PETSc work vectors tmpx and tmpy point to the arrays in the PVODE vectors
    163:   */
    164:   VecPlaceArray(tmpx,N_VGetData(y));
    165:   if (ierr) {
    166:     (*PetscErrorPrintf)("TSFunction_PVode:Could not place array. Error code %d",(int)ierr);
    167:   }
    168:   VecPlaceArray(tmpy,N_VGetData(ydot));
    169:   if (ierr) {
    170:     (*PetscErrorPrintf)("TSFunction_PVode:Could not place array. Error code %d",(int)ierr);
    171:   }

    173:   /* now compute the right hand side function */
    174:   TSComputeRHSFunction(ts,t,tmpx,tmpy);
    175:   if (ierr) {
    176:     (*PetscErrorPrintf)("TSFunction_PVode:Could not compute RHS function. Error code %d",(int)ierr);
    177:   }
    178: }

    180: /*
    181:        TSStep_PVode_Nonlinear - Calls PVode to integrate the ODE.

    183:     Contributed by: Liyang Xu
    184: */
    187: /*
    188:     TSStep_PVode_Nonlinear -
    189:
    190:    steps - number of time steps
    191:    time - time that integrater is  terminated.

    193: */
    194: PetscErrorCode TSStep_PVode_Nonlinear(TS ts,int *steps,double *time)
    195: {
    196:   TS_PVode  *cvode = (TS_PVode*)ts->data;
    197:   Vec       sol = ts->vec_sol;
    198:   int       ierr;
    199:   int i,max_steps = ts->max_steps,flag;
    200:   double    t,tout;
    201:   realtype  *tmp;

    204:   /* initialize the number of steps */
    205:   *steps = -ts->steps;
    206:   TSMonitor(ts,ts->steps,ts->ptime,sol);

    208:   /* call CVSpgmr to use GMRES as the linear solver. */
    209:   /* setup the ode integrator with the given preconditioner */
    210:   CVSpgmr(cvode->mem,LEFT,cvode->gtype,cvode->restart,cvode->linear_tol,TSPrecond_PVode,TSPSolve_PVode,ts,0,0);

    212:   tout = ts->max_time;
    213:   for (i=0; i<max_steps; i++) {
    214:     if (ts->ptime >= tout) break;
    215:     VecGetArray(ts->vec_sol,&tmp);
    216:     N_VSetData(tmp,cvode->y);
    217:     flag = CVode(cvode->mem,tout,cvode->y,&t,ONE_STEP);
    218:     cvode->nonlinear_solves += cvode->iopt[NNI];
    219:     VecRestoreArray(ts->vec_sol,PETSC_NULL);
    220:     if (flag != SUCCESS) SETERRQ(PETSC_ERR_LIB,"PVODE failed");

    222:     if (t > tout && cvode->exact_final_time) {
    223:       /* interpolate to final requested time */
    224:       VecGetArray(ts->vec_sol,&tmp);
    225:       N_VSetData(tmp,cvode->y);
    226:       flag = CVodeDky(cvode->mem,tout,0,cvode->y);
    227:       VecRestoreArray(ts->vec_sol,PETSC_NULL);
    228:       if (flag != SUCCESS) SETERRQ(PETSC_ERR_LIB,"PVODE interpolation to final time failed");
    229:       t = tout;
    230:     }

    232:     ts->time_step = t - ts->ptime;
    233:     ts->ptime     = t;

    235:     /*
    236:        copy the solution from cvode->y to cvode->update and sol
    237:     */
    238:     VecPlaceArray(cvode->w1,N_VGetData(cvode->y));
    239:     VecCopy(cvode->w1,cvode->update);
    240:     VecCopy(cvode->update,sol);
    241:
    242:     ts->steps++;
    243:     TSMonitor(ts,ts->steps,t,sol);
    244:     ts->nonlinear_its = cvode->iopt[NNI];
    245:     ts->linear_its    = cvode->iopt[SPGMR_NLI];
    246:   }

    248:   *steps           += ts->steps;
    249:   *time             = t;

    251:   return(0);
    252: }

    254: /*

    256:     Contributed by: Liyang Xu
    257: */
    260: PetscErrorCode TSDestroy_PVode(TS ts)
    261: {
    262:   TS_PVode *cvode = (TS_PVode*)ts->data;

    266:   if (cvode->pmat)   {MatDestroy(cvode->pmat);}
    267:   if (cvode->pc)     {PCDestroy(cvode->pc);}
    268:   if (cvode->update) {VecDestroy(cvode->update);}
    269:   if (cvode->func)   {VecDestroy(cvode->func);}
    270:   if (cvode->rhs)    {VecDestroy(cvode->rhs);}
    271:   if (cvode->w1)     {VecDestroy(cvode->w1);}
    272:   if (cvode->w2)     {VecDestroy(cvode->w2);}
    273:   MPI_Comm_free(&(cvode->comm_pvode));
    274:   PetscFree(cvode);
    275:   return(0);
    276: }

    278: /*

    280:     Contributed by: Liyang Xu
    281: */
    284: PetscErrorCode TSSetUp_PVode_Nonlinear(TS ts)
    285: {
    286:   TS_PVode    *cvode = (TS_PVode*)ts->data;
    288:   int M,locsize;
    289:   M_Env       machEnv;
    290:   realtype    *tmp;

    293:   PCSetFromOptions(cvode->pc);
    294:   /* get the vector size */
    295:   VecGetSize(ts->vec_sol,&M);
    296:   VecGetLocalSize(ts->vec_sol,&locsize);

    298:   /* allocate the memory for machEnv */
    299:   /* machEnv = PVInitMPI(cvode->>comm_pvode,locsize,M);  */
    300:   machEnv = M_EnvInit_Parallel(cvode->comm_pvode, locsize, M, 0, 0);

    303:   /* allocate the memory for N_Vec y */
    304:   cvode->y         = N_VNew(M,machEnv);
    305:   VecGetArray(ts->vec_sol,&tmp);
    306:   N_VSetData(tmp,cvode->y);
    307:   VecRestoreArray(ts->vec_sol,PETSC_NULL);

    309:   /* initializing vector update and func */
    310:   VecDuplicate(ts->vec_sol,&cvode->update);
    311:   VecDuplicate(ts->vec_sol,&cvode->func);
    312:   PetscLogObjectParent(ts,cvode->update);
    313:   PetscLogObjectParent(ts,cvode->func);

    315:   /*
    316:       Create work vectors for the TSPSolve_PVode() routine. Note these are
    317:     allocated with zero space arrays because the actual array space is provided
    318:     by PVode and set using VecPlaceArray().
    319:   */
    320:   VecCreateMPIWithArray(ts->comm,locsize,PETSC_DECIDE,0,&cvode->w1);
    321:   VecCreateMPIWithArray(ts->comm,locsize,PETSC_DECIDE,0,&cvode->w2);
    322:   PetscLogObjectParent(ts,cvode->w1);
    323:   PetscLogObjectParent(ts,cvode->w2);

    325:   /* allocate memory for PVode */
    326:   VecGetArray(ts->vec_sol,&tmp);
    327:   N_VSetData(tmp,cvode->y);
    328:   cvode->mem = CVodeMalloc(M,TSFunction_PVode,ts->ptime,cvode->y,cvode->cvode_type,
    329:                            NEWTON,SS,&cvode->reltol,&cvode->abstol,ts,NULL,TRUE,cvode->iopt,
    330:                            cvode->ropt,machEnv);
    331:   VecRestoreArray(ts->vec_sol,PETSC_NULL);
    332:   return(0);
    333: }

    335: /*

    337:     Contributed by: Liyang Xu
    338: */
    341: PetscErrorCode TSSetFromOptions_PVode_Nonlinear(TS ts)
    342: {
    343:   TS_PVode   *cvode = (TS_PVode*)ts->data;
    345:   int indx;
    346:   const char *btype[] = {"bdf","adams"},*otype[] = {"modified","unmodified"};
    347:   PetscTruth flag;

    350:   PetscOptionsHead("PVODE ODE solver options");
    351:     PetscOptionsEList("-ts_pvode_type","Scheme","TSPVodeSetType",btype,2,"bdf",&indx,&flag);
    352:     if (flag) {
    353:       TSPVodeSetType(ts,(TSPVodeType)indx);
    354:     }
    355:     PetscOptionsEList("-ts_pvode_gramschmidt_type","Type of orthogonalization","TSPVodeSetGramSchmidtType",otype,2,"unmodified",&indx,&flag);
    356:     if (flag) {
    357:       TSPVodeSetGramSchmidtType(ts,(TSPVodeGramSchmidtType)indx);
    358:     }
    359:     PetscOptionsReal("-ts_pvode_atol","Absolute tolerance for convergence","TSPVodeSetTolerance",cvode->abstol,&cvode->abstol,PETSC_NULL);
    360:     PetscOptionsReal("-ts_pvode_rtol","Relative tolerance for convergence","TSPVodeSetTolerance",cvode->reltol,&cvode->reltol,PETSC_NULL);
    361:     PetscOptionsReal("-ts_pvode_linear_tolerance","Convergence tolerance for linear solve","TSPVodeSetLinearTolerance",cvode->linear_tol,&cvode->linear_tol,&flag);
    362:     PetscOptionsInt("-ts_pvode_gmres_restart","Number of GMRES orthogonalization directions","TSPVodeSetGMRESRestart",cvode->restart,&cvode->restart,&flag);
    363:     PetscOptionsName("-ts_pvode_not_exact_final_time","Allow PVODE to stop near the final time, not exactly on it","TSPVodeSetExactFinalTime",&cvode->exact_final_time);
    364:   PetscOptionsTail();

    366:   return(0);
    367: }

    369: /*

    371:     Contributed by: Liyang Xu
    372: */
    375: PetscErrorCode TSPrintHelp_PVode(TS ts,char *p)
    376: {

    380:   (*PetscHelpPrintf)(ts->comm," Options for TSPVODE integrater:\n");
    381:   (*PetscHelpPrintf)(ts->comm," -ts_pvode_type <bdf,adams>: integration approach\n",p);
    382:   (*PetscHelpPrintf)(ts->comm," -ts_pvode_atol aabs: absolute tolerance of ODE solution\n",p);
    383:   (*PetscHelpPrintf)(ts->comm," -ts_pvode_rtol rel: relative tolerance of ODE solution\n",p);
    384:   (*PetscHelpPrintf)(ts->comm," -ts_pvode_gramschmidt_type <unmodified,modified>\n");
    385:   (*PetscHelpPrintf)(ts->comm," -ts_pvode_gmres_restart <restart_size> (also max. GMRES its)\n");
    386:   (*PetscHelpPrintf)(ts->comm," -ts_pvode_linear_tolerance <tol>\n");
    387:   (*PetscHelpPrintf)(ts->comm," -ts_pvode_not_exact_final_time\n");

    389:   return(0);
    390: }

    392: /*

    394:     Contributed by: Liyang Xu
    395: */
    398: PetscErrorCode TSView_PVode(TS ts,PetscViewer viewer)
    399: {
    400:   TS_PVode   *cvode = (TS_PVode*)ts->data;
    402:   char       *type;
    403:   PetscTruth iascii,isstring;

    406:   if (cvode->cvode_type == PVODE_ADAMS) {type = "Adams";}
    407:   else                                  {type = "BDF: backward differentiation formula";}

    409:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    410:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
    411:   if (iascii) {
    412:     PetscViewerASCIIPrintf(viewer,"PVode integrater does not use SNES!\n");
    413:     PetscViewerASCIIPrintf(viewer,"PVode integrater type %s\n",type);
    414:     PetscViewerASCIIPrintf(viewer,"PVode abs tol %g rel tol %g\n",cvode->abstol,cvode->reltol);
    415:     PetscViewerASCIIPrintf(viewer,"PVode linear solver tolerance factor %g\n",cvode->linear_tol);
    416:     PetscViewerASCIIPrintf(viewer,"PVode GMRES max iterations (same as restart in PVODE) %D\n",cvode->restart);
    417:     if (cvode->gtype == PVODE_MODIFIED_GS) {
    418:       PetscViewerASCIIPrintf(viewer,"PVode using modified Gram-Schmidt for orthogonalization in GMRES\n");
    419:     } else {
    420:       PetscViewerASCIIPrintf(viewer,"PVode using unmodified (classical) Gram-Schmidt for orthogonalization in GMRES\n");
    421:     }
    422:   } else if (isstring) {
    423:     PetscViewerStringSPrintf(viewer,"Pvode type %s",type);
    424:   } else {
    425:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by TS PVode",((PetscObject)viewer)->type_name);
    426:   }
    427:   PetscViewerASCIIPushTab(viewer);
    428:   PCView(cvode->pc,viewer);
    429:   PetscViewerASCIIPopTab(viewer);

    431:   return(0);
    432: }

    435: /* --------------------------------------------------------------------------*/
    439: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeSetType_PVode(TS ts,TSPVodeType type)
    440: {
    441:   TS_PVode *cvode = (TS_PVode*)ts->data;
    442:
    444:   cvode->cvode_type = type;
    445:   return(0);
    446: }

    452: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeSetGMRESRestart_PVode(TS ts,int restart)
    453: {
    454:   TS_PVode *cvode = (TS_PVode*)ts->data;
    455:
    457:   cvode->restart = restart;
    458:   return(0);
    459: }

    465: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeSetLinearTolerance_PVode(TS ts,double tol)
    466: {
    467:   TS_PVode *cvode = (TS_PVode*)ts->data;
    468:
    470:   cvode->linear_tol = tol;
    471:   return(0);
    472: }

    478: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeSetGramSchmidtType_PVode(TS ts,TSPVodeGramSchmidtType type)
    479: {
    480:   TS_PVode *cvode = (TS_PVode*)ts->data;
    481:
    483:   cvode->gtype = type;

    485:   return(0);
    486: }

    492: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeSetTolerance_PVode(TS ts,double aabs,double rel)
    493: {
    494:   TS_PVode *cvode = (TS_PVode*)ts->data;
    495:
    497:   if (aabs != PETSC_DECIDE) cvode->abstol = aabs;
    498:   if (rel != PETSC_DECIDE)  cvode->reltol = rel;
    499:   return(0);
    500: }

    506: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeGetPC_PVode(TS ts,PC *pc)
    507: {
    508:   TS_PVode *cvode = (TS_PVode*)ts->data;

    511:   *pc = cvode->pc;

    513:   return(0);
    514: }

    520: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeGetIterations_PVode(TS ts,int *nonlin,int *lin)
    521: {
    522:   TS_PVode *cvode = (TS_PVode*)ts->data;
    523:
    525:   if (nonlin) *nonlin = cvode->nonlinear_solves;
    526:   if (lin)    *lin    = cvode->linear_solves;
    527:   return(0);
    528: }
    530:
    534: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeSetExactFinalTime_PVode(TS ts,PetscTruth s)
    535: {
    536:   TS_PVode *cvode = (TS_PVode*)ts->data;
    537:
    539:   cvode->exact_final_time = s;
    540:   return(0);
    541: }
    543: /* -------------------------------------------------------------------------------------------*/

    547: /*@C
    548:    TSPVodeGetIterations - Gets the number of nonlinear and linear iterations used so far by PVode.

    550:    Not Collective

    552:    Input parameters:
    553: .    ts     - the time-step context

    555:    Output Parameters:
    556: +   nonlin - number of nonlinear iterations
    557: -   lin    - number of linear iterations

    559:    Level: advanced

    561:    Notes:
    562:     These return the number since the creation of the TS object

    564: .keywords: non-linear iterations, linear iterations

    566: .seealso: TSPVodeSetType(), TSPVodeSetGMRESRestart(),
    567:           TSPVodeSetLinearTolerance(), TSPVodeSetGramSchmidtType(), TSPVodeSetTolerance(),
    568:           TSPVodeGetIterations(), TSPVodeSetType(), TSPVodeSetGMRESRestart(),
    569:           TSPVodeSetLinearTolerance(), TSPVodeSetTolerance(), TSPVodeGetPC(),
    570:           TSPVodeSetExactFinalTime()

    572: @*/
    573: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeGetIterations(TS ts,int *nonlin,int *lin)
    574: {
    575:   PetscErrorCode ierr,(*f)(TS,int*,int*);
    576:
    578:   PetscObjectQueryFunction((PetscObject)ts,"TSPVodeGetIterations_C",(void (**)(void))&f);
    579:   if (f) {
    580:     (*f)(ts,nonlin,lin);
    581:   }
    582:   return(0);
    583: }

    587: /*@
    588:    TSPVodeSetType - Sets the method that PVode will use for integration.

    590:    Collective on TS

    592:    Input parameters:
    593: +    ts     - the time-step context
    594: -    type   - one of  PVODE_ADAMS or PVODE_BDF

    596:     Contributed by: Liyang Xu

    598:    Level: intermediate

    600: .keywords: Adams, backward differentiation formula

    602: .seealso: TSPVodeGetIterations(),  TSPVodeSetGMRESRestart(),
    603:           TSPVodeSetLinearTolerance(), TSPVodeSetGramSchmidtType(), TSPVodeSetTolerance(),
    604:           TSPVodeGetIterations(), TSPVodeSetType(), TSPVodeSetGMRESRestart(),
    605:           TSPVodeSetLinearTolerance(), TSPVodeSetTolerance(), TSPVodeGetPC(),
    606:           TSPVodeSetExactFinalTime()
    607: @*/
    608: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeSetType(TS ts,TSPVodeType type)
    609: {
    610:   PetscErrorCode ierr,(*f)(TS,TSPVodeType);
    611:
    613:   PetscObjectQueryFunction((PetscObject)ts,"TSPVodeSetType_C",(void (**)(void))&f);
    614:   if (f) {
    615:     (*f)(ts,type);
    616:   }
    617:   return(0);
    618: }

    622: /*@
    623:    TSPVodeSetGMRESRestart - Sets the dimension of the Krylov space used by
    624:        GMRES in the linear solver in PVODE. PVODE DOES NOT use restarted GMRES so
    625:        this is ALSO the maximum number of GMRES steps that will be used.

    627:    Collective on TS

    629:    Input parameters:
    630: +    ts      - the time-step context
    631: -    restart - number of direction vectors (the restart size).

    633:    Level: advanced

    635: .keywords: GMRES, restart

    637: .seealso: TSPVodeGetIterations(), TSPVodeSetType(),
    638:           TSPVodeSetLinearTolerance(), TSPVodeSetGramSchmidtType(), TSPVodeSetTolerance(),
    639:           TSPVodeGetIterations(), TSPVodeSetType(), TSPVodeSetGMRESRestart(),
    640:           TSPVodeSetLinearTolerance(), TSPVodeSetTolerance(), TSPVodeGetPC(),
    641:           TSPVodeSetExactFinalTime()

    643: @*/
    644: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeSetGMRESRestart(TS ts,int restart)
    645: {
    646:   PetscErrorCode ierr,(*f)(TS,int);

    649:   PetscObjectQueryFunction((PetscObject)ts,"TSPVodeSetGMRESRestart_C",(void (**)(void))&f);
    650:   if (f) {
    651:     (*f)(ts,restart);
    652:   }

    654:   return(0);
    655: }

    659: /*@
    660:    TSPVodeSetLinearTolerance - Sets the tolerance used to solve the linear
    661:        system by PVODE.

    663:    Collective on TS

    665:    Input parameters:
    666: +    ts     - the time-step context
    667: -    tol    - the factor by which the tolerance on the nonlinear solver is
    668:              multiplied to get the tolerance on the linear solver, .05 by default.

    670:    Level: advanced

    672: .keywords: GMRES, linear convergence tolerance, PVODE

    674: .seealso: TSPVodeGetIterations(), TSPVodeSetType(), TSPVodeSetGMRESRestart(),
    675:           TSPVodeSetGramSchmidtType(), TSPVodeSetTolerance(),
    676:           TSPVodeGetIterations(), TSPVodeSetType(), TSPVodeSetGMRESRestart(),
    677:           TSPVodeSetLinearTolerance(), TSPVodeSetTolerance(), TSPVodeGetPC(),
    678:           TSPVodeSetExactFinalTime()

    680: @*/
    681: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeSetLinearTolerance(TS ts,double tol)
    682: {
    683:   PetscErrorCode ierr,(*f)(TS,double);
    684:
    686:   PetscObjectQueryFunction((PetscObject)ts,"TSPVodeSetLinearTolerance_C",(void (**)(void))&f);
    687:   if (f) {
    688:     (*f)(ts,tol);
    689:   }
    690:   return(0);
    691: }

    695: /*@
    696:    TSPVodeSetGramSchmidtType - Sets type of orthogonalization used
    697:         in GMRES method by PVODE linear solver.

    699:    Collective on TS

    701:    Input parameters:
    702: +    ts  - the time-step context
    703: -    type - either PVODE_MODIFIED_GS or PVODE_CLASSICAL_GS

    705:    Level: advanced

    707: .keywords: PVode, orthogonalization

    709: .seealso: TSPVodeGetIterations(), TSPVodeSetType(), TSPVodeSetGMRESRestart(),
    710:           TSPVodeSetLinearTolerance(),  TSPVodeSetTolerance(),
    711:           TSPVodeGetIterations(), TSPVodeSetType(), TSPVodeSetGMRESRestart(),
    712:           TSPVodeSetLinearTolerance(), TSPVodeSetTolerance(), TSPVodeGetPC(),
    713:           TSPVodeSetExactFinalTime()

    715: @*/
    716: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeSetGramSchmidtType(TS ts,TSPVodeGramSchmidtType type)
    717: {
    718:   PetscErrorCode ierr,(*f)(TS,TSPVodeGramSchmidtType);
    719:
    721:   PetscObjectQueryFunction((PetscObject)ts,"TSPVodeSetGramSchmidtType_C",(void (**)(void))&f);
    722:   if (f) {
    723:     (*f)(ts,type);
    724:   }
    725:   return(0);
    726: }

    730: /*@
    731:    TSPVodeSetTolerance - Sets the absolute and relative tolerance used by
    732:                          PVode for error control.

    734:    Collective on TS

    736:    Input parameters:
    737: +    ts  - the time-step context
    738: .    aabs - the absolute tolerance
    739: -    rel - the relative tolerance

    741:     Contributed by: Liyang Xu

    743:      See the Cvode/Pvode users manual for exact details on these parameters. Essentially
    744:     these regulate the size of the error for a SINGLE timestep.

    746:    Level: intermediate

    748: .keywords: PVode, tolerance

    750: .seealso: TSPVodeGetIterations(), TSPVodeSetType(), TSPVodeSetGMRESRestart(),
    751:           TSPVodeSetLinearTolerance(), TSPVodeSetGramSchmidtType(),
    752:           TSPVodeGetIterations(), TSPVodeSetType(), TSPVodeSetGMRESRestart(),
    753:           TSPVodeSetLinearTolerance(), TSPVodeSetTolerance(), TSPVodeGetPC(),
    754:           TSPVodeSetExactFinalTime()

    756: @*/
    757: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeSetTolerance(TS ts,double aabs,double rel)
    758: {
    759:   PetscErrorCode ierr,(*f)(TS,double,double);
    760:
    762:   PetscObjectQueryFunction((PetscObject)ts,"TSPVodeSetTolerance_C",(void (**)(void))&f);
    763:   if (f) {
    764:     (*f)(ts,aabs,rel);
    765:   }
    766:   return(0);
    767: }

    771: /*@
    772:    TSPVodeGetPC - Extract the PC context from a time-step context for PVode.

    774:    Input Parameter:
    775: .    ts - the time-step context

    777:    Output Parameter:
    778: .    pc - the preconditioner context

    780:    Level: advanced

    782:     Contributed by: Liyang Xu

    784: .seealso: TSPVodeGetIterations(), TSPVodeSetType(), TSPVodeSetGMRESRestart(),
    785:           TSPVodeSetLinearTolerance(), TSPVodeSetGramSchmidtType(), TSPVodeSetTolerance(),
    786:           TSPVodeGetIterations(), TSPVodeSetType(), TSPVodeSetGMRESRestart(),
    787:           TSPVodeSetLinearTolerance(), TSPVodeSetTolerance()
    788: @*/
    789: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeGetPC(TS ts,PC *pc)
    790: {
    791:   PetscErrorCode ierr,(*f)(TS,PC *);

    794:   PetscObjectQueryFunction((PetscObject)ts,"TSPVodeGetPC_C",(void (**)(void))&f);
    795:   if (f) {
    796:     (*f)(ts,pc);
    797:   } else {
    798:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"TS must be of PVode type to extract the PC");
    799:   }

    801:   return(0);
    802: }

    806: /*@
    807:    TSPVodeSetExactFinalTime - Determines if PVode interpolates solution to the
    808:       exact final time requested by the user or just returns it at the final time
    809:       it computed. (Defaults to true).

    811:    Input Parameter:
    812: +   ts - the time-step context
    813: -   ft - PETSC_TRUE if interpolates, else PETSC_FALSE

    815:    Level: beginner

    817: .seealso:TSPVodeGetIterations(), TSPVodeSetType(), TSPVodeSetGMRESRestart(),
    818:           TSPVodeSetLinearTolerance(), TSPVodeSetGramSchmidtType(), TSPVodeSetTolerance(),
    819:           TSPVodeGetIterations(), TSPVodeSetType(), TSPVodeSetGMRESRestart(),
    820:           TSPVodeSetLinearTolerance(), TSPVodeSetTolerance(), TSPVodeGetPC()
    821: @*/
    822: PetscErrorCode PETSCTS_DLLEXPORT TSPVodeSetExactFinalTime(TS ts,PetscTruth ft)
    823: {
    824:   PetscErrorCode ierr,(*f)(TS,PetscTruth);

    827:   PetscObjectQueryFunction((PetscObject)ts,"TSPVodeSetExactFinalTime_C",(void (**)(void))&f);
    828:   if (f) {
    829:     (*f)(ts,ft);
    830:   }

    832:   return(0);
    833: }

    835: /* -------------------------------------------------------------------------------------------*/
    836: /*MC
    837:       TS_PVode - ODE solver using the LLNL CVODE/PVODE package (now called SUNDIALS)

    839:    Options Database:
    840: +    -ts_pvode_type <bdf,adams>
    841: .    -ts_pvode_gramschmidt_type <modified, classical> - type of orthogonalization inside GMRES
    842: .    -ts_pvode_atol <tol> - Absolute tolerance for convergence
    843: .    -ts_pvode_rtol <tol> - Relative tolerance for convergence
    844: .    -ts_pvode_linear_tolerance <tol>
    845: .    -ts_pvode_gmres_restart <restart> - Number of GMRES orthogonalization directions
    846: -    -ts_pvode_not_exact_final_time -Allow PVODE to stop near the final time, not exactly on it

    848:     Notes: This uses its own nonlinear solver and Krylov method so PETSc SNES and KSP options do not apply
    849:            only PETSc PC options

    851:     Contributed by: Liyang Xu

    853:     Level: beginner

    855: .seealso:  TSCreate(), TS, TSSetType(), TSPVodeSetType(), TSPVodeSetGMRESRestart(), TSPVodeSetLinearTolerance(),
    856:            TSPVodeSetGramSchmidtType(), TSPVodeSetTolerance(), TSPVodeGetPC(), TSPVodeGetIterations(), TSPVodeSetExactFinalTime()

    858: M*/
    862: PetscErrorCode PETSCTS_DLLEXPORT TSCreate_PVode(TS ts)
    863: {
    864:   TS_PVode *cvode;

    868:   ts->ops->destroy         = TSDestroy_PVode;
    869:   ts->ops->view            = TSView_PVode;

    871:   if (ts->problem_type != TS_NONLINEAR) {
    872:     SETERRQ(PETSC_ERR_SUP,"Only support for nonlinear problems");
    873:   }
    874:   ts->ops->setup           = TSSetUp_PVode_Nonlinear;
    875:   ts->ops->step            = TSStep_PVode_Nonlinear;
    876:   ts->ops->setfromoptions  = TSSetFromOptions_PVode_Nonlinear;

    878:   PetscNew(TS_PVode,&cvode);
    879:   PCCreate(ts->comm,&cvode->pc);
    880:   PetscLogObjectParent(ts,cvode->pc);
    881:   ts->data          = (void*)cvode;
    882:   cvode->cvode_type = BDF;
    883:   cvode->gtype      = PVODE_UNMODIFIED_GS;
    884:   cvode->restart    = 5;
    885:   cvode->linear_tol = .05;

    887:   cvode->exact_final_time = PETSC_TRUE;

    889:   MPI_Comm_dup(ts->comm,&(cvode->comm_pvode));
    890:   /* set tolerance for PVode */
    891:   cvode->abstol = 1e-6;
    892:   cvode->reltol = 1e-6;

    894:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPVodeSetType_C","TSPVodeSetType_PVode",
    895:                     TSPVodeSetType_PVode);
    896:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPVodeSetGMRESRestart_C",
    897:                     "TSPVodeSetGMRESRestart_PVode",
    898:                     TSPVodeSetGMRESRestart_PVode);
    899:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPVodeSetLinearTolerance_C",
    900:                     "TSPVodeSetLinearTolerance_PVode",
    901:                      TSPVodeSetLinearTolerance_PVode);
    902:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPVodeSetGramSchmidtType_C",
    903:                     "TSPVodeSetGramSchmidtType_PVode",
    904:                      TSPVodeSetGramSchmidtType_PVode);
    905:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPVodeSetTolerance_C",
    906:                     "TSPVodeSetTolerance_PVode",
    907:                      TSPVodeSetTolerance_PVode);
    908:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPVodeGetPC_C",
    909:                     "TSPVodeGetPC_PVode",
    910:                      TSPVodeGetPC_PVode);
    911:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPVodeGetIterations_C",
    912:                     "TSPVodeGetIterations_PVode",
    913:                      TSPVodeGetIterations_PVode);
    914:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPVodeSetExactFinalTime_C",
    915:                     "TSPVodeSetExactFinalTime_PVode",
    916:                      TSPVodeSetExactFinalTime_PVode);
    917:   return(0);
    918: }