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

      1: static char help[] = "Solves DAE with integrator only on non-algebraic terms \n";

      3: #include <petscts.h>

      5: /*
      6:         \dot{U} = f(U,V)
      7:         F(U,V)  = 0

      9:     Same as ex6.c and ex7.c except calls the ARKIMEX integrator on the entire DAE
     10: */

     15: /*
     16:    f(U,V) = U + V

     18: */
     19: PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F)
     20: {

     24:   VecWAXPY(F,1.0,U,V);
     25:   return(0);
     26: }

     30: /*
     31:    F(U,V) = U - V

     33: */
     34: PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F)
     35: {

     39:   VecWAXPY(F,-1.0,V,U);
     40:   return(0);
     41: }

     44: typedef struct {
     45:   Vec            U,V;
     46:   Vec            UF,VF;
     47:   VecScatter     scatterU,scatterV;
     48:   PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec);
     49:   PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec);
     50: } AppCtx;

     52: extern PetscErrorCode TSFunctionRHS(TS,PetscReal,Vec,Vec,void*);
     53: extern PetscErrorCode TSFunctionI(TS,PetscReal,Vec,Vec,Vec,void*);

     57: int main(int argc,char **argv)
     58: {
     60:   AppCtx         ctx;
     61:   TS             ts;
     62:   Vec            tsrhs,UV;
     63:   IS             is;
     64:   PetscInt       I;
     65:   PetscMPIInt    rank;

     68:   PetscInitialize(&argc,&argv,(char*)0,help);
     69:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
     70:   TSCreate(PETSC_COMM_WORLD,&ts);
     71:   TSSetProblemType(ts,TS_NONLINEAR);
     72:   TSSetType(ts,TSROSW);
     73:   TSSetFromOptions(ts);
     74:   VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&tsrhs);
     75:   VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&UV);
     76:   TSSetRHSFunction(ts,tsrhs,TSFunctionRHS,&ctx);
     77:   TSSetIFunction(ts,NULL,TSFunctionI,&ctx);
     78:   ctx.f = f;
     79:   ctx.F = F;

     81:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.U);
     82:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.V);
     83:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.UF);
     84:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.VF);
     85:   I    = 2*rank;
     86:   ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is);
     87:   VecScatterCreate(ctx.U,NULL,UV,is,&ctx.scatterU);
     88:   ISDestroy(&is);
     89:   I    = 2*rank + 1;
     90:   ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is);
     91:   VecScatterCreate(ctx.V,NULL,UV,is,&ctx.scatterV);
     92:   ISDestroy(&is);

     94:   VecSet(UV,1.0);
     95:   TSSolve(ts,UV);
     96:   VecDestroy(&tsrhs);
     97:   VecDestroy(&UV);
     98:   VecDestroy(&ctx.U);
     99:   VecDestroy(&ctx.V);
    100:   VecDestroy(&ctx.UF);
    101:   VecDestroy(&ctx.VF);
    102:   VecScatterDestroy(&ctx.scatterU);
    103:   VecScatterDestroy(&ctx.scatterV);
    104:   TSDestroy(&ts);
    105:   PetscFinalize();
    106:   return 0;
    107: }

    111: /*
    112:    Defines the RHS function that is passed to the time-integrator.

    114: */
    115: PetscErrorCode TSFunctionRHS(TS ts,PetscReal t,Vec UV,Vec F,void *actx)
    116: {
    117:   AppCtx         *ctx = (AppCtx*)actx;

    121:   VecSet(F,0.0);
    122:   VecScatterBegin(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE);
    123:   VecScatterEnd(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE);
    124:   VecScatterBegin(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE);
    125:   VecScatterEnd(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE);
    126:   (*ctx->f)(t,ctx->U,ctx->V,ctx->UF);
    127:   VecScatterBegin(ctx->scatterU,ctx->UF,F,INSERT_VALUES,SCATTER_FORWARD);
    128:   VecScatterEnd(ctx->scatterU,ctx->UF,F,INSERT_VALUES,SCATTER_FORWARD);
    129:   return(0);
    130: }

    134: /*
    135:    Defines the nonlinear function that is passed to the time-integrator

    137: */
    138: PetscErrorCode TSFunctionI(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx)
    139: {
    140:   AppCtx         *ctx = (AppCtx*)actx;

    144:   VecCopy(UVdot,F);
    145:   VecScatterBegin(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE);
    146:   VecScatterEnd(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE);
    147:   VecScatterBegin(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE);
    148:   VecScatterEnd(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE);
    149:   (*ctx->F)(t,ctx->U,ctx->V,ctx->VF);
    150:   VecScatterBegin(ctx->scatterV,ctx->VF,F,INSERT_VALUES,SCATTER_FORWARD);
    151:   VecScatterEnd(ctx->scatterV,ctx->VF,F,INSERT_VALUES,SCATTER_FORWARD);
    152:   return(0);
    153: }