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


 11: /*
 12:    f(U,V) = U + V
 13: */
 14: PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F)
 15: {

 19:   VecWAXPY(F,1.0,U,V);
 20:   return(0);
 21: }

 23: /*
 24:    F(U,V) = U - V
 25: */
 26: PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F)
 27: {

 31:   VecWAXPY(F,-1.0,V,U);
 32:   return(0);
 33: }

 35: typedef struct {
 36:   PetscReal      t;
 37:   SNES           snes;
 38:   Vec            U,V;
 39:   PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec);
 40:   PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec);
 41: } AppCtx;

 43: extern PetscErrorCode TSFunction(TS,PetscReal,Vec,Vec,void*);
 44: extern PetscErrorCode SNESFunction(SNES,Vec,Vec,void*);

 46: int main(int argc,char **argv)
 47: {
 49:   AppCtx         ctx;
 50:   TS             ts;
 51:   Vec            tsrhs,U;

 53:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 54:   TSCreate(PETSC_COMM_WORLD,&ts);
 55:   TSSetProblemType(ts,TS_NONLINEAR);
 56:   TSSetType(ts,TSEULER);
 57:   TSSetFromOptions(ts);
 58:   VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&tsrhs);
 59:   VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&U);
 60:   TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx);
 61:   TSSetMaxTime(ts,1.0);
 62:   ctx.f = f;

 64:   SNESCreate(PETSC_COMM_WORLD,&ctx.snes);
 65:   SNESSetFromOptions(ctx.snes);
 66:   SNESSetFunction(ctx.snes,NULL,SNESFunction,&ctx);
 67:   SNESSetJacobian(ctx.snes,NULL,NULL,SNESComputeJacobianDefault,&ctx);
 68:   ctx.F = F;
 69:   VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&ctx.V);

 71:   VecSet(U,1.0);
 72:   TSSolve(ts,U);

 74:   VecDestroy(&ctx.V);
 75:   VecDestroy(&tsrhs);
 76:   VecDestroy(&U);
 77:   SNESDestroy(&ctx.snes);
 78:   TSDestroy(&ts);
 79:   PetscFinalize();
 80:   return ierr;
 81: }

 83: /*
 84:    Defines the RHS function that is passed to the time-integrator.

 86:    Solves F(U,V) for V and then computes f(U,V)
 87: */
 88: PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
 89: {
 90:   AppCtx         *ctx = (AppCtx*)actx;

 94:   ctx->t = t;
 95:   ctx->U = U;
 96:   SNESSolve(ctx->snes,NULL,ctx->V);
 97:   (*ctx->f)(t,U,ctx->V,F);
 98:   return(0);
 99: }

101: /*
102:    Defines the nonlinear function that is passed to the nonlinear solver
103: */
104: PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx)
105: {
106:   AppCtx         *ctx = (AppCtx*)actx;

110:   (*ctx->F)(ctx->t,ctx->U,V,F);
111:   return(0);
112: }

114: /*TEST

116:     test:
117:       args:  -ts_monitor -ts_view

119: TEST*/