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 TSROSW integrator on the entire DAE
 10: */


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

 16: */
 17: PetscErrorCode f(PetscReal t,Vec UV,Vec F)
 18: {
 19:   PetscErrorCode    ierr;
 20:   const PetscScalar *u,*v;
 21:   PetscScalar       *f;
 22:   PetscInt          n,i;

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

 36: /*
 37:    F(U,V) = U - V

 39: */
 40: PetscErrorCode F(PetscReal t,Vec UV,Vec F)
 41: {
 42:   PetscErrorCode    ierr;
 43:   const PetscScalar *u,*v;
 44:   PetscScalar       *f;
 45:   PetscInt          n,i;

 48:   VecGetLocalSize(UV,&n);
 49:   n    = n/2;
 50:   VecGetArrayRead(UV,&u);
 51:   v    = u + n;
 52:   VecGetArrayWrite(F,&f);
 53:   f    = f + n;
 54:   for (i=0; i<n; i++) f[i] = u[i] - v[i];
 55:   f    = f - n;
 56:   VecRestoreArrayRead(UV,&u);
 57:   VecRestoreArrayWrite(F,&f);
 58:   return(0);
 59: }

 61: typedef struct {
 62:   PetscErrorCode (*f)(PetscReal,Vec,Vec);
 63:   PetscErrorCode (*F)(PetscReal,Vec,Vec);
 64: } AppCtx;

 66: extern PetscErrorCode TSFunctionRHS(TS,PetscReal,Vec,Vec,void*);
 67: extern PetscErrorCode TSFunctionI(TS,PetscReal,Vec,Vec,Vec,void*);

 69: int main(int argc,char **argv)
 70: {
 72:   AppCtx         ctx;
 73:   TS             ts;
 74:   Vec            tsrhs,UV;

 76:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 77:   TSCreate(PETSC_COMM_WORLD,&ts);
 78:   TSSetProblemType(ts,TS_NONLINEAR);
 79:   TSSetType(ts,TSROSW);
 80:   TSSetFromOptions(ts);
 81:   VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&tsrhs);
 82:   VecDuplicate(tsrhs,&UV);
 83:   TSSetRHSFunction(ts,tsrhs,TSFunctionRHS,&ctx);
 84:   TSSetIFunction(ts,NULL,TSFunctionI,&ctx);
 85:   TSSetMaxTime(ts,1.0);
 86:   ctx.f = f;
 87:   ctx.F = F;

 89:   VecSet(UV,1.0);
 90:   TSSolve(ts,UV);
 91:   VecDestroy(&tsrhs);
 92:   VecDestroy(&UV);
 93:   TSDestroy(&ts);
 94:   PetscFinalize();
 95:   return ierr;
 96: }

 98: /*
 99:    Defines the RHS function that is passed to the time-integrator.
100: */
101: PetscErrorCode TSFunctionRHS(TS ts,PetscReal t,Vec UV,Vec F,void *actx)
102: {
103:   AppCtx         *ctx = (AppCtx*)actx;

107:   VecSet(F,0.0);
108:   (*ctx->f)(t,UV,F);
109:   return(0);
110: }

112: /*
113:    Defines the nonlinear function that is passed to the time-integrator
114: */
115: PetscErrorCode TSFunctionI(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx)
116: {
117:   AppCtx         *ctx = (AppCtx*)actx;

121:   VecCopy(UVdot,F);
122:   (*ctx->F)(t,UV,F);
123:   return(0);
124: }

126: /*TEST

128:     test:
129:       args:  -ts_view

131:     test:
132:       suffix: 2
133:       args: -snes_lag_jacobian 2 -ts_view

135: TEST*/