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

  9:     Same as ex6.c and ex7.c except calls the TSROSW integrator on the entire DAE
 10: */

 12: /*
 13:    f(U,V) = U + V

 15: */
 16: PetscErrorCode f(PetscReal t, Vec UV, Vec F)
 17: {
 18:   const PetscScalar *u, *v;
 19:   PetscScalar       *f;
 20:   PetscInt           n, i;

 22:   PetscFunctionBeginUser;
 23:   PetscCall(VecGetLocalSize(UV, &n));
 24:   n = n / 2;
 25:   PetscCall(VecGetArrayRead(UV, &u));
 26:   v = u + n;
 27:   PetscCall(VecGetArrayWrite(F, &f));
 28:   for (i = 0; i < n; i++) f[i] = u[i] + v[i];
 29:   PetscCall(VecRestoreArrayRead(UV, &u));
 30:   PetscCall(VecRestoreArrayWrite(F, &f));
 31:   PetscFunctionReturn(PETSC_SUCCESS);
 32: }

 34: /*
 35:    F(U,V) = U - V

 37: */
 38: PetscErrorCode F(PetscReal t, Vec UV, Vec F)
 39: {
 40:   const PetscScalar *u, *v;
 41:   PetscScalar       *f;
 42:   PetscInt           n, i;

 44:   PetscFunctionBeginUser;
 45:   PetscCall(VecGetLocalSize(UV, &n));
 46:   n = n / 2;
 47:   PetscCall(VecGetArrayRead(UV, &u));
 48:   v = u + n;
 49:   PetscCall(VecGetArrayWrite(F, &f));
 50:   f = f + n;
 51:   for (i = 0; i < n; i++) f[i] = u[i] - v[i];
 52:   f = f - n;
 53:   PetscCall(VecRestoreArrayRead(UV, &u));
 54:   PetscCall(VecRestoreArrayWrite(F, &f));
 55:   PetscFunctionReturn(PETSC_SUCCESS);
 56: }

 58: typedef struct {
 59:   PetscErrorCode (*f)(PetscReal, Vec, Vec);
 60:   PetscErrorCode (*F)(PetscReal, Vec, Vec);
 61: } AppCtx;

 63: extern PetscErrorCode TSFunctionRHS(TS, PetscReal, Vec, Vec, void *);
 64: extern PetscErrorCode TSFunctionI(TS, PetscReal, Vec, Vec, Vec, void *);

 66: int main(int argc, char **argv)
 67: {
 68:   AppCtx ctx;
 69:   TS     ts;
 70:   Vec    tsrhs, UV;

 72:   PetscFunctionBeginUser;
 73:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 74:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
 75:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
 76:   PetscCall(TSSetType(ts, TSROSW));
 77:   PetscCall(TSSetFromOptions(ts));
 78:   PetscCall(VecCreateMPI(PETSC_COMM_WORLD, 2, PETSC_DETERMINE, &tsrhs));
 79:   PetscCall(VecDuplicate(tsrhs, &UV));
 80:   PetscCall(TSSetRHSFunction(ts, tsrhs, TSFunctionRHS, &ctx));
 81:   PetscCall(TSSetIFunction(ts, NULL, TSFunctionI, &ctx));
 82:   PetscCall(TSSetMaxTime(ts, 1.0));
 83:   ctx.f = f;
 84:   ctx.F = F;

 86:   PetscCall(VecSet(UV, 1.0));
 87:   PetscCall(TSSolve(ts, UV));
 88:   PetscCall(VecDestroy(&tsrhs));
 89:   PetscCall(VecDestroy(&UV));
 90:   PetscCall(TSDestroy(&ts));
 91:   PetscCall(PetscFinalize());
 92:   return 0;
 93: }

 95: /*
 96:    Defines the RHS function that is passed to the time-integrator.
 97: */
 98: PetscErrorCode TSFunctionRHS(TS ts, PetscReal t, Vec UV, Vec F, void *actx)
 99: {
100:   AppCtx *ctx = (AppCtx *)actx;

102:   PetscFunctionBeginUser;
103:   PetscCall(VecSet(F, 0.0));
104:   PetscCall((*ctx->f)(t, UV, F));
105:   PetscFunctionReturn(PETSC_SUCCESS);
106: }

108: /*
109:    Defines the nonlinear function that is passed to the time-integrator
110: */
111: PetscErrorCode TSFunctionI(TS ts, PetscReal t, Vec UV, Vec UVdot, Vec F, void *actx)
112: {
113:   AppCtx *ctx = (AppCtx *)actx;

115:   PetscFunctionBeginUser;
116:   PetscCall(VecCopy(UVdot, F));
117:   PetscCall((*ctx->F)(t, UV, F));
118:   PetscFunctionReturn(PETSC_SUCCESS);
119: }

121: /*TEST

123:     test:
124:       args: -ts_view

126:     test:
127:       suffix: 2
128:       args: -snes_lag_jacobian 2 -ts_view

130: TEST*/