Actual source code: ex6.c

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

 10: /*
 11:    f(U,V) = U + V
 12: */
 13: PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F)
 14: {
 16:   VecWAXPY(F,1.0,U,V);
 17:   return 0;
 18: }

 20: /*
 21:    F(U,V) = U - V
 22: */
 23: PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F)
 24: {
 26:   VecWAXPY(F,-1.0,V,U);
 27:   return 0;
 28: }

 30: typedef struct {
 31:   PetscReal      t;
 32:   SNES           snes;
 33:   Vec            U,V;
 34:   PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec);
 35:   PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec);
 36: } AppCtx;

 38: extern PetscErrorCode TSFunction(TS,PetscReal,Vec,Vec,void*);
 39: extern PetscErrorCode SNESFunction(SNES,Vec,Vec,void*);

 41: int main(int argc,char **argv)
 42: {
 43:   AppCtx         ctx;
 44:   TS             ts;
 45:   Vec            tsrhs,U;

 47:   PetscInitialize(&argc,&argv,(char*)0,help);
 48:   TSCreate(PETSC_COMM_WORLD,&ts);
 49:   TSSetProblemType(ts,TS_NONLINEAR);
 50:   TSSetType(ts,TSEULER);
 51:   TSSetFromOptions(ts);
 52:   VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&tsrhs);
 53:   VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&U);
 54:   TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx);
 55:   TSSetMaxTime(ts,1.0);
 56:   ctx.f = f;

 58:   SNESCreate(PETSC_COMM_WORLD,&ctx.snes);
 59:   SNESSetFromOptions(ctx.snes);
 60:   SNESSetFunction(ctx.snes,NULL,SNESFunction,&ctx);
 61:   SNESSetJacobian(ctx.snes,NULL,NULL,SNESComputeJacobianDefault,&ctx);
 62:   ctx.F = F;
 63:   VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&ctx.V);

 65:   VecSet(U,1.0);
 66:   TSSolve(ts,U);

 68:   VecDestroy(&ctx.V);
 69:   VecDestroy(&tsrhs);
 70:   VecDestroy(&U);
 71:   SNESDestroy(&ctx.snes);
 72:   TSDestroy(&ts);
 73:   PetscFinalize();
 74:   return 0;
 75: }

 77: /*
 78:    Defines the RHS function that is passed to the time-integrator.

 80:    Solves F(U,V) for V and then computes f(U,V)
 81: */
 82: PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
 83: {
 84:   AppCtx         *ctx = (AppCtx*)actx;

 87:   ctx->t = t;
 88:   ctx->U = U;
 89:   SNESSolve(ctx->snes,NULL,ctx->V);
 90:   (*ctx->f)(t,U,ctx->V,F);
 91:   return 0;
 92: }

 94: /*
 95:    Defines the nonlinear function that is passed to the nonlinear solver
 96: */
 97: PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx)
 98: {
 99:   AppCtx         *ctx = (AppCtx*)actx;

102:   (*ctx->F)(ctx->t,ctx->U,V,F);
103:   return 0;
104: }

106: /*TEST

108:     test:
109:       args:  -ts_monitor -ts_view

111: TEST*/