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


 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: }

 28: /*
 29:    F(U,V) = U - V

 31: */
 32: PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F)
 33: {

 37:   VecWAXPY(F,-1.0,V,U);
 38:   return(0);
 39: }

 41: typedef struct {
 42:   PetscReal      t;
 43:   SNES           snes;
 44:   Vec            U,V;
 45:   PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec);
 46:   PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec);
 47: } AppCtx;

 49: extern PetscErrorCode TSFunction(TS,PetscReal,Vec,Vec,void*);
 50: extern PetscErrorCode SNESFunction(SNES,Vec,Vec,void*);

 54: int main(int argc,char **argv)
 55: {
 57:   AppCtx         ctx;
 58:   TS             ts;
 59:   Vec            tsrhs,U;

 61:   PetscInitialize(&argc,&argv,(char*)0,help);
 62:   TSCreate(PETSC_COMM_WORLD,&ts);
 63:   TSSetProblemType(ts,TS_NONLINEAR);
 64:   TSSetType(ts,TSEULER);
 65:   TSSetFromOptions(ts);
 66:   VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&tsrhs);
 67:   VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&U);
 68:   TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx);
 69:   ctx.f = f;

 71:   SNESCreate(PETSC_COMM_WORLD,&ctx.snes);
 72:   SNESSetFromOptions(ctx.snes);
 73:   SNESSetFunction(ctx.snes,NULL,SNESFunction,&ctx);
 74:   SNESSetJacobian(ctx.snes,NULL,NULL,SNESComputeJacobianDefault,&ctx);
 75:   ctx.F = F;
 76:   VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&ctx.V);

 78:   VecSet(U,1.0);
 79:   TSSolve(ts,U);

 81:   VecDestroy(&ctx.V);
 82:   VecDestroy(&tsrhs);
 83:   VecDestroy(&U);
 84:   SNESDestroy(&ctx.snes);
 85:   TSDestroy(&ts);
 86:   PetscFinalize();
 87:   return 0;
 88: }

 92: /*
 93:    Defines the RHS function that is passed to the time-integrator. 

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

 97: */
 98: PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
 99: {
100:   AppCtx         *ctx = (AppCtx*)actx;

104:   ctx->t = t;
105:   ctx->U = U;
106:   SNESSolve(ctx->snes,NULL,ctx->V);
107:   (*ctx->f)(t,U,ctx->V,F);
108:   return(0);
109: }

113: /*
114:    Defines the nonlinear function that is passed to the nonlinear solver

116: */
117: PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx)
118: {
119:   AppCtx         *ctx = (AppCtx*)actx;

123:   (*ctx->F)(ctx->t,ctx->U,V,F);
124:   return(0);
125: }