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: */


 13: /*
 14:    f(U,V) = U + V

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

 22:   VecWAXPY(F,1.0,U,V);
 23:   return(0);
 24: }

 26: /*
 27:    F(U,V) = U - V

 29: */
 30: PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F)
 31: {

 35:   VecWAXPY(F,-1.0,V,U);
 36:   return(0);
 37: }


 40: typedef struct {
 41:   Vec            U,V;
 42:   Vec            UF,VF;
 43:   VecScatter     scatterU,scatterV;
 44:   PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec);
 45:   PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec);
 46: } AppCtx;

 48: extern PetscErrorCode TSFunctionRHS(TS,PetscReal,Vec,Vec,void*);
 49: extern PetscErrorCode TSFunctionI(TS,PetscReal,Vec,Vec,Vec,void*);

 51: int main(int argc,char **argv)
 52: {
 54:   AppCtx         ctx;
 55:   TS             ts;
 56:   Vec            tsrhs,UV;
 57:   IS             is;
 58:   PetscInt       I;
 59:   PetscMPIInt    rank;


 62:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 63:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 64:   TSCreate(PETSC_COMM_WORLD,&ts);
 65:   TSSetProblemType(ts,TS_NONLINEAR);
 66:   TSSetType(ts,TSROSW);
 67:   TSSetFromOptions(ts);
 68:   VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&tsrhs);
 69:   VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&UV);
 70:   TSSetRHSFunction(ts,tsrhs,TSFunctionRHS,&ctx);
 71:   TSSetIFunction(ts,NULL,TSFunctionI,&ctx);
 72:   ctx.f = f;
 73:   ctx.F = F;

 75:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.U);
 76:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.V);
 77:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.UF);
 78:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.VF);
 79:   I    = 2*rank;
 80:   ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is);
 81:   VecScatterCreate(ctx.U,NULL,UV,is,&ctx.scatterU);
 82:   ISDestroy(&is);
 83:   I    = 2*rank + 1;
 84:   ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is);
 85:   VecScatterCreate(ctx.V,NULL,UV,is,&ctx.scatterV);
 86:   ISDestroy(&is);

 88:   VecSet(UV,1.0);
 89:   TSSolve(ts,UV);
 90:   VecDestroy(&tsrhs);
 91:   VecDestroy(&UV);
 92:   VecDestroy(&ctx.U);
 93:   VecDestroy(&ctx.V);
 94:   VecDestroy(&ctx.UF);
 95:   VecDestroy(&ctx.VF);
 96:   VecScatterDestroy(&ctx.scatterU);
 97:   VecScatterDestroy(&ctx.scatterV);
 98:   TSDestroy(&ts);
 99:   PetscFinalize();
100:   return ierr;
101: }

103: /*
104:    Defines the RHS function that is passed to the time-integrator. 

106: */
107: PetscErrorCode TSFunctionRHS(TS ts,PetscReal t,Vec UV,Vec F,void *actx)
108: {
109:   AppCtx         *ctx = (AppCtx*)actx;

113:   VecSet(F,0.0);
114:   VecScatterBegin(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE);
115:   VecScatterEnd(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE);
116:   VecScatterBegin(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE);
117:   VecScatterEnd(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE);
118:   (*ctx->f)(t,ctx->U,ctx->V,ctx->UF);
119:   VecScatterBegin(ctx->scatterU,ctx->UF,F,INSERT_VALUES,SCATTER_FORWARD);
120:   VecScatterEnd(ctx->scatterU,ctx->UF,F,INSERT_VALUES,SCATTER_FORWARD);
121:   return(0);
122: }

124: /*
125:    Defines the nonlinear function that is passed to the time-integrator

127: */
128: PetscErrorCode TSFunctionI(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx)
129: {
130:   AppCtx         *ctx = (AppCtx*)actx;

134:   VecCopy(UVdot,F);
135:   VecScatterBegin(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE);
136:   VecScatterEnd(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE);
137:   VecScatterBegin(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE);
138:   VecScatterEnd(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE);
139:   (*ctx->F)(t,ctx->U,ctx->V,ctx->VF);
140:   VecScatterBegin(ctx->scatterV,ctx->VF,F,INSERT_VALUES,SCATTER_FORWARD);
141:   VecScatterEnd(ctx->scatterV,ctx->VF,F,INSERT_VALUES,SCATTER_FORWARD);
142:   return(0);
143: }