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: {
 15:   PetscFunctionBeginUser;
 16:   PetscCall(VecWAXPY(F, 1.0, U, V));
 17:   PetscFunctionReturn(PETSC_SUCCESS);
 18: }

 20: /*
 21:    F(U,V) = U - V
 22: */
 23: PetscErrorCode F(PetscReal t, Vec U, Vec V, Vec F)
 24: {
 25:   PetscFunctionBeginUser;
 26:   PetscCall(VecWAXPY(F, -1.0, V, U));
 27:   PetscFunctionReturn(PETSC_SUCCESS);
 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:   PetscFunctionBeginUser;
 48:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 49:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
 50:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
 51:   PetscCall(TSSetType(ts, TSEULER));
 52:   PetscCall(TSSetFromOptions(ts));
 53:   PetscCall(VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, 1, &tsrhs));
 54:   PetscCall(VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, 1, &U));
 55:   PetscCall(TSSetRHSFunction(ts, tsrhs, TSFunction, &ctx));
 56:   PetscCall(TSSetMaxTime(ts, 1.0));
 57:   ctx.f = f;

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

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

 69:   PetscCall(VecDestroy(&ctx.V));
 70:   PetscCall(VecDestroy(&tsrhs));
 71:   PetscCall(VecDestroy(&U));
 72:   PetscCall(SNESDestroy(&ctx.snes));
 73:   PetscCall(TSDestroy(&ts));
 74:   PetscCall(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;

 87:   PetscFunctionBeginUser;
 88:   ctx->t = t;
 89:   ctx->U = U;
 90:   PetscCall(SNESSolve(ctx->snes, NULL, ctx->V));
 91:   PetscCall((*ctx->f)(t, U, ctx->V, F));
 92:   PetscFunctionReturn(PETSC_SUCCESS);
 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;

102:   PetscFunctionBeginUser;
103:   PetscCall((*ctx->F)(ctx->t, ctx->U, V, F));
104:   PetscFunctionReturn(PETSC_SUCCESS);
105: }

107: /*TEST

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

112: TEST*/