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 UV,Vec F)
 20: {
 21:   PetscErrorCode    ierr;
 22:   const PetscScalar *u,*v;
 23:   PetscScalar       *f;
 24:   PetscInt          n,i;

 27:   VecGetLocalSize(UV,&n);
 28:   n    = n/2;
 29:   VecGetArrayRead(UV,&u);
 30:   v    = u + n;
 31:   VecGetArray(F,&f);
 32:   for (i=0; i<n; i++) f[i] = u[i] + v[i];
 33:   VecRestoreArrayRead(UV,&u);
 34:   VecRestoreArray(F,&f);
 35:   return(0);
 36: }

 40: /*
 41:    F(U,V) = U - V

 43: */
 44: PetscErrorCode F(PetscReal t,Vec UV,Vec F)
 45: {
 46:   PetscErrorCode    ierr;
 47:   const PetscScalar *u,*v;
 48:   PetscScalar       *f;
 49:   PetscInt          n,i;

 52:   VecGetLocalSize(UV,&n);
 53:   n    = n/2;
 54:   VecGetArrayRead(UV,&u);
 55:   v    = u + n;
 56:   VecGetArray(F,&f);
 57:   f    = f + n;
 58:   for (i=0; i<n; i++) f[i] = u[i] - v[i];
 59:   f    = f - n;
 60:   VecRestoreArrayRead(UV,&u);
 61:   VecRestoreArray(F,&f);
 62:   return(0);
 63: }

 65: typedef struct {
 66:   PetscErrorCode (*f)(PetscReal,Vec,Vec);
 67:   PetscErrorCode (*F)(PetscReal,Vec,Vec);
 68: } AppCtx;

 70: extern PetscErrorCode TSFunctionRHS(TS,PetscReal,Vec,Vec,void*);
 71: extern PetscErrorCode TSFunctionI(TS,PetscReal,Vec,Vec,Vec,void*);

 75: int main(int argc,char **argv)
 76: {
 78:   AppCtx         ctx;
 79:   TS             ts;
 80:   Vec            tsrhs,UV;

 82:   PetscInitialize(&argc,&argv,(char*)0,help);
 83:   TSCreate(PETSC_COMM_WORLD,&ts);
 84:   TSSetProblemType(ts,TS_NONLINEAR);
 85:   TSSetType(ts,TSROSW);
 86:   TSSetFromOptions(ts);
 87:   VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&tsrhs);
 88:   VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&UV);
 89:   TSSetRHSFunction(ts,tsrhs,TSFunctionRHS,&ctx);
 90:   TSSetIFunction(ts,NULL,TSFunctionI,&ctx);
 91:   ctx.f = f;
 92:   ctx.F = F;

 94:   VecSet(UV,1.0);
 95:   TSSolve(ts,UV);
 96:   VecDestroy(&tsrhs);
 97:   VecDestroy(&UV);
 98:   TSDestroy(&ts);
 99:   PetscFinalize();
100:   return 0;
101: }

105: /*
106:    Defines the RHS function that is passed to the time-integrator. 



110: */
111: PetscErrorCode TSFunctionRHS(TS ts,PetscReal t,Vec UV,Vec F,void *actx)
112: {
113:   AppCtx         *ctx = (AppCtx*)actx;

117:   VecSet(F,0.0);
118:   (*ctx->f)(t,UV,F);
119:   return(0);
120: }

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:   (*ctx->F)(t,UV,F);
136:   return(0);
137: }