Actual source code: ex10.c
petsc-3.4.2 2013-07-02

      1: static char help[] = "Simple wrapper object to solve DAE of the form:\n\
      2:                              \\dot{U} = f(U,V)\n\
      3:                              F(U,V) = 0\n\n";

      5: #include <petscts.h>

      7: /* ----------------------------------------------------------------------------*/

      9: typedef struct _p_TSDAESimple *TSDAESimple;
     10: struct _p_TSDAESimple {
     11:   MPI_Comm       comm;
     12:   PetscErrorCode (*setfromoptions)(TSDAESimple);
     13:   PetscErrorCode (*solve)(TSDAESimple,Vec);
     14:   PetscErrorCode (*destroy)(TSDAESimple);
     15:   Vec            U,V;
     16:   PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec,void*);
     17:   PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec,void*);
     18:   void           *fctx,*Fctx;
     19:   void           *data;
     20: };

     24: PetscErrorCode TSDAESimpleCreate(MPI_Comm comm,TSDAESimple *tsdae)
     25: {

     29:   PetscNew(struct _p_TSDAESimple,tsdae);
     30:   (*tsdae)->comm = comm;
     31:   return(0);
     32: }

     36: PetscErrorCode TSDAESimpleSetRHSFunction(TSDAESimple tsdae,Vec U,PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
     37: {

     41:   tsdae->f    = f;
     42:   tsdae->U    = U;
     43:   PetscObjectReference((PetscObject)U);
     44:   tsdae->fctx = ctx;
     45:   return(0);
     46: }

     50: PetscErrorCode TSDAESimpleSetIFunction(TSDAESimple tsdae,Vec V,PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
     51: {

     55:   tsdae->F    = F;
     56:   tsdae->V    = V;
     57:   PetscObjectReference((PetscObject)V);
     58:   tsdae->Fctx = ctx;
     59:   return(0);
     60: }

     64: PetscErrorCode TSDAESimpleDestroy(TSDAESimple *tsdae)
     65: {

     69:   (*(*tsdae)->destroy)(*tsdae);
     70:   VecDestroy(&(*tsdae)->U);
     71:   VecDestroy(&(*tsdae)->V);
     72:   PetscFree(*tsdae);
     73:   return(0);
     74: }

     78: PetscErrorCode TSDAESimpleSolve(TSDAESimple tsdae,Vec Usolution)
     79: {

     83:   (*tsdae->solve)(tsdae,Usolution);
     84:   return(0);
     85: }

     89: PetscErrorCode TSDAESimpleSetFromOptions(TSDAESimple tsdae)
     90: {

     94:   (*tsdae->setfromoptions)(tsdae);
     95:   return(0);
     96: }

     98: /* ----------------------------------------------------------------------------*/
     99: /*
    100:       Integrates the system by integrating the reduced ODE system and solving the
    101:    algebraic constraints at each stage with a separate SNES solve.
    102: */

    104: typedef struct {
    105:   PetscReal t;
    106:   TS        ts;
    107:   SNES      snes;
    108:   Vec       U;
    109: } TSDAESimple_Reduced;

    113: /*
    114:    Defines the RHS function that is passed to the time-integrator.

    116:    Solves F(U,V) for V and then computes f(U,V)

    118: */
    119: PetscErrorCode TSDAESimple_Reduced_TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
    120: {
    121:   TSDAESimple         tsdae = (TSDAESimple)actx;
    122:   TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
    123:   PetscErrorCode      ierr;

    126:   red->t = t;
    127:   red->U = U;
    128:   SNESSolve(red->snes,NULL,tsdae->V);
    129:   (*tsdae->f)(t,U,tsdae->V,F,tsdae->fctx);
    130:   return(0);
    131: }

    135: /*
    136:    Defines the nonlinear function that is passed to the nonlinear solver

    138: */
    139: PetscErrorCode TSDAESimple_Reduced_SNESFunction(SNES snes,Vec V,Vec F,void *actx)
    140: {
    141:   TSDAESimple         tsdae = (TSDAESimple)actx;
    142:   TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
    143:   PetscErrorCode      ierr;

    146:   (*tsdae->F)(red->t,red->U,V,F,tsdae->Fctx);
    147:   return(0);
    148: }

    153: PetscErrorCode TSDAESimpleSolve_Reduced(TSDAESimple tsdae,Vec U)
    154: {
    155:   PetscErrorCode      ierr;
    156:   TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;

    159:   TSSolve(red->ts,U);
    160:   return(0);
    161: }

    165: PetscErrorCode TSDAESimpleSetFromOptions_Reduced(TSDAESimple tsdae)
    166: {
    167:   PetscErrorCode      ierr;
    168:   TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;

    171:   TSSetFromOptions(red->ts);
    172:   SNESSetFromOptions(red->snes);
    173:   return(0);
    174: }

    178: PetscErrorCode TSDAESimpleDestroy_Reduced(TSDAESimple tsdae)
    179: {
    180:   PetscErrorCode      ierr;
    181:   TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;

    184:   TSDestroy(&red->ts);
    185:   SNESDestroy(&red->snes);
    186:   PetscFree(red);
    187:   return(0);
    188: }

    192: PetscErrorCode TSDAESimpleSetUp_Reduced(TSDAESimple tsdae)
    193: {
    194:   PetscErrorCode      ierr;
    195:   TSDAESimple_Reduced *red;
    196:   Vec                 tsrhs;

    199:   PetscMalloc(sizeof(TSDAESimple_Reduced),&red);
    200:   tsdae->data = red;

    202:   tsdae->setfromoptions = TSDAESimpleSetFromOptions_Reduced;
    203:   tsdae->solve          = TSDAESimpleSolve_Reduced;
    204:   tsdae->destroy        = TSDAESimpleDestroy_Reduced;

    206:   TSCreate(tsdae->comm,&red->ts);
    207:   TSSetProblemType(red->ts,TS_NONLINEAR);
    208:   TSSetType(red->ts,TSEULER);
    209:   VecDuplicate(tsdae->U,&tsrhs);
    210:   TSSetRHSFunction(red->ts,tsrhs,TSDAESimple_Reduced_TSFunction,tsdae);
    211:   VecDestroy(&tsrhs);

    213:   SNESCreate(tsdae->comm,&red->snes);
    214:   SNESSetOptionsPrefix(red->snes,"tsdaesimple_");
    215:   SNESSetFunction(red->snes,NULL,TSDAESimple_Reduced_SNESFunction,tsdae);
    216:   SNESSetJacobian(red->snes,NULL,NULL,SNESComputeJacobianDefault,tsdae);
    217:   return(0);
    218: }

    221: /* ----------------------------------------------------------------------------*/

    223: /*
    224:       Integrates the system by integrating directly the entire DAE system
    225: */

    227: typedef struct {
    228:   TS         ts;
    229:   Vec        UV,UF,VF;
    230:   VecScatter scatterU,scatterV;
    231: } TSDAESimple_Full;

    235: /*
    236:    Defines the RHS function that is passed to the time-integrator.

    238:    f(U,V)
    239:    0

    241: */
    242: PetscErrorCode TSDAESimple_Full_TSRHSFunction(TS ts,PetscReal t,Vec UV,Vec F,void *actx)
    243: {
    244:   TSDAESimple      tsdae = (TSDAESimple)actx;
    245:   TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
    246:   PetscErrorCode   ierr;

    249:   VecSet(F,0.0);
    250:   VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
    251:   VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
    252:   VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
    253:   VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
    254:   (*tsdae->f)(t,tsdae->U,tsdae->V,full->UF,tsdae->fctx);
    255:   VecScatterBegin(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);
    256:   VecScatterEnd(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);
    257:   return(0);
    258: }

    262: /*
    263:    Defines the nonlinear function that is passed to the nonlinear solver

    265:    \dot{U}
    266:    F(U,V)

    268: */
    269: PetscErrorCode TSDAESimple_Full_TSIFunction(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx)
    270: {
    271:   TSDAESimple       tsdae = (TSDAESimple)actx;
    272:   TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
    273:   PetscErrorCode    ierr;

    276:   VecCopy(UVdot,F);
    277:   VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
    278:   VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
    279:   VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
    280:   VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
    281:   (*tsdae->F)(t,tsdae->U,tsdae->V,full->VF,tsdae->Fctx);
    282:   VecScatterBegin(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);
    283:   VecScatterEnd(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);
    284:   return(0);
    285: }

    290: PetscErrorCode TSDAESimpleSolve_Full(TSDAESimple tsdae,Vec U)
    291: {
    292:   PetscErrorCode   ierr;
    293:   TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;

    296:   VecSet(full->UV,1.0);
    297:   VecScatterBegin(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);
    298:   VecScatterEnd(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);
    299:   TSSolve(full->ts,full->UV);
    300:   VecScatterBegin(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);
    301:   VecScatterEnd(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);
    302:   return(0);
    303: }

    307: PetscErrorCode TSDAESimpleSetFromOptions_Full(TSDAESimple tsdae)
    308: {
    309:   PetscErrorCode   ierr;
    310:   TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;

    313:   TSSetFromOptions(full->ts);
    314:   return(0);
    315: }

    319: PetscErrorCode TSDAESimpleDestroy_Full(TSDAESimple tsdae)
    320: {
    321:   PetscErrorCode   ierr;
    322:   TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;

    325:   TSDestroy(&full->ts);
    326:   VecDestroy(&full->UV);
    327:   VecDestroy(&full->UF);
    328:   VecDestroy(&full->VF);
    329:   VecScatterDestroy(&full->scatterU);
    330:   VecScatterDestroy(&full->scatterV);
    331:   PetscFree(full);
    332:   return(0);
    333: }

    337: PetscErrorCode TSDAESimpleSetUp_Full(TSDAESimple tsdae)
    338: {
    339:   PetscErrorCode   ierr;
    340:   TSDAESimple_Full *full;
    341:   Vec              tsrhs;
    342:   PetscInt         nU,nV,UVstart;
    343:   IS               is;

    346:   PetscMalloc(sizeof(TSDAESimple_Full),&full);
    347:   tsdae->data = full;

    349:   tsdae->setfromoptions = TSDAESimpleSetFromOptions_Full;
    350:   tsdae->solve          = TSDAESimpleSolve_Full;
    351:   tsdae->destroy        = TSDAESimpleDestroy_Full;

    353:   TSCreate(tsdae->comm,&full->ts);
    354:   TSSetProblemType(full->ts,TS_NONLINEAR);
    355:   TSSetType(full->ts,TSROSW);
    356:   VecDuplicate(tsdae->U,&full->UF);
    357:   VecDuplicate(tsdae->V,&full->VF);

    359:   VecGetLocalSize(tsdae->U,&nU);
    360:   VecGetLocalSize(tsdae->V,&nV);
    361:   VecCreateMPI(tsdae->comm,nU+nV,PETSC_DETERMINE,&tsrhs);
    362:   VecDuplicate(tsrhs,&full->UV);

    364:   VecGetOwnershipRange(tsrhs,&UVstart,NULL);
    365:   ISCreateStride(tsdae->comm,nU,UVstart,1,&is);
    366:   VecScatterCreate(tsdae->U,NULL,tsrhs,is,&full->scatterU);
    367:   ISDestroy(&is);
    368:   ISCreateStride(tsdae->comm,nV,UVstart+nU,1,&is);
    369:   VecScatterCreate(tsdae->V,NULL,tsrhs,is,&full->scatterV);
    370:   ISDestroy(&is);

    372:   TSSetRHSFunction(full->ts,tsrhs,TSDAESimple_Full_TSRHSFunction,tsdae);
    373:   TSSetIFunction(full->ts,NULL,TSDAESimple_Full_TSIFunction,tsdae);
    374:   VecDestroy(&tsrhs);
    375:   return(0);
    376: }

    379: /* ----------------------------------------------------------------------------*/

    384: /*
    385:    Simple example:   f(U,V) = U + V

    387: */
    388: PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
    389: {

    393:   VecWAXPY(F,1.0,U,V);
    394:   return(0);
    395: }

    399: /*
    400:    Simple example: F(U,V) = U - V

    402: */
    403: PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
    404: {

    408:   VecWAXPY(F,-1.0,V,U);
    409:   return(0);
    410: }

    414: int main(int argc,char **argv)
    415: {
    417:   TSDAESimple    tsdae;
    418:   Vec            U,V,Usolution;

    420:   PetscInitialize(&argc,&argv,(char*)0,help);
    421:   TSDAESimpleCreate(PETSC_COMM_WORLD,&tsdae);

    423:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&U);
    424:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&V);
    425:   TSDAESimpleSetRHSFunction(tsdae,U,f,NULL);
    426:   TSDAESimpleSetIFunction(tsdae,V,F,NULL);

    428:   VecDuplicate(U,&Usolution);
    429:   VecSet(Usolution,1.0);

    431:   /*  TSDAESimpleSetUp_Full(tsdae); */
    432:   TSDAESimpleSetUp_Reduced(tsdae);

    434:   TSDAESimpleSetFromOptions(tsdae);
    435:   TSDAESimpleSolve(tsdae,Usolution);
    436:   TSDAESimpleDestroy(&tsdae);

    438:   VecDestroy(&U);
    439:   VecDestroy(&Usolution);
    440:   VecDestroy(&V);
    441:   PetscFinalize();
    442:   return 0;
    443: }