Actual source code: ex8.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

  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;

 23:   VecGetLocalSize(UV, &n);
 24:   n = n / 2;
 25:   VecGetArrayRead(UV, &u);
 26:   v = u + n;
 27:   VecGetArrayWrite(F, &f);
 28:   for (i = 0; i < n; i++) f[i] = u[i] + v[i];
 29:   VecRestoreArrayRead(UV, &u);
 30:   VecRestoreArrayWrite(F, &f);
 31:   return 0;
 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;

 45:   VecGetLocalSize(UV, &n);
 46:   n = n / 2;
 47:   VecGetArrayRead(UV, &u);
 48:   v = u + n;
 49:   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:   VecRestoreArrayRead(UV, &u);
 54:   VecRestoreArrayWrite(F, &f);
 55:   return 0;
 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;

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

 86:   VecSet(UV, 1.0);
 87:   TSSolve(ts, UV);
 88:   VecDestroy(&tsrhs);
 89:   VecDestroy(&UV);
 90:   TSDestroy(&ts);
 91:   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;

103:   VecSet(F, 0.0);
104:   (*ctx->f)(t, UV, F);
105:   return 0;
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;

116:   VecCopy(UVdot, F);
117:   (*ctx->F)(t, UV, F);
118:   return 0;
119: }

121: /*TEST

123:     test:
124:       args:  -ts_view

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

130: TEST*/