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;

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

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

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

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

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

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

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

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

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

107: /*TEST

109:     test:
110:       args:  -ts_monitor -ts_view

112: TEST*/