Actual source code: ex7.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 except the user provided functions take input values as a single vector instead of two vectors
 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
 36: */
 37: PetscErrorCode F(PetscReal t, Vec UV, Vec F)
 38: {
 39:   const PetscScalar *u, *v;
 40:   PetscScalar       *f;
 41:   PetscInt           n, i;

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

 55: typedef struct {
 56:   PetscReal  t;
 57:   SNES       snes;
 58:   Vec        UV, V;
 59:   VecScatter scatterU, scatterV;
 60:   PetscErrorCode (*f)(PetscReal, Vec, Vec);
 61:   PetscErrorCode (*F)(PetscReal, Vec, Vec);
 62: } AppCtx;

 64: extern PetscErrorCode TSFunction(TS, PetscReal, Vec, Vec, void *);
 65: extern PetscErrorCode SNESFunction(SNES, Vec, Vec, void *);

 67: int main(int argc, char **argv)
 68: {
 69:   AppCtx      ctx;
 70:   TS          ts;
 71:   Vec         tsrhs, U;
 72:   IS          is;
 73:   PetscInt    i;
 74:   PetscMPIInt rank;

 77:   PetscInitialize(&argc, &argv, (char *)0, help);
 78:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 79:   TSCreate(PETSC_COMM_WORLD, &ts);
 80:   TSSetProblemType(ts, TS_NONLINEAR);
 81:   TSSetType(ts, TSEULER);
 82:   TSSetFromOptions(ts);
 83:   VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &tsrhs);
 84:   VecDuplicate(tsrhs, &U);
 85:   TSSetRHSFunction(ts, tsrhs, TSFunction, &ctx);
 86:   TSSetMaxTime(ts, 1.0);
 87:   ctx.f = f;

 89:   SNESCreate(PETSC_COMM_WORLD, &ctx.snes);
 90:   SNESSetFromOptions(ctx.snes);
 91:   SNESSetFunction(ctx.snes, NULL, SNESFunction, &ctx);
 92:   SNESSetJacobian(ctx.snes, NULL, NULL, SNESComputeJacobianDefault, &ctx);
 93:   ctx.F = F;
 94:   VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.V);

 96:   /* Create scatters to move between separate U and V representation and UV representation of solution */
 97:   VecCreateMPI(PETSC_COMM_WORLD, 2, PETSC_DETERMINE, &ctx.UV);
 98:   i = 2 * rank;
 99:   ISCreateGeneral(PETSC_COMM_WORLD, 1, &i, PETSC_COPY_VALUES, &is);
100:   VecScatterCreate(U, NULL, ctx.UV, is, &ctx.scatterU);
101:   ISDestroy(&is);
102:   i = 2 * rank + 1;
103:   ISCreateGeneral(PETSC_COMM_WORLD, 1, &i, PETSC_COPY_VALUES, &is);
104:   VecScatterCreate(ctx.V, NULL, ctx.UV, is, &ctx.scatterV);
105:   ISDestroy(&is);

107:   VecSet(U, 1.0);
108:   TSSolve(ts, U);

110:   VecDestroy(&ctx.V);
111:   VecDestroy(&ctx.UV);
112:   VecScatterDestroy(&ctx.scatterU);
113:   VecScatterDestroy(&ctx.scatterV);
114:   VecDestroy(&tsrhs);
115:   VecDestroy(&U);
116:   SNESDestroy(&ctx.snes);
117:   TSDestroy(&ts);
118:   PetscFinalize();
119:   return 0;
120: }

122: /*
123:    Defines the RHS function that is passed to the time-integrator.

125:    Solves F(U,V) for V and then computes f(U,V)
126: */
127: PetscErrorCode TSFunction(TS ts, PetscReal t, Vec U, Vec F, void *actx)
128: {
129:   AppCtx *ctx = (AppCtx *)actx;

132:   ctx->t = t;
133:   VecScatterBegin(ctx->scatterU, U, ctx->UV, INSERT_VALUES, SCATTER_FORWARD);
134:   VecScatterEnd(ctx->scatterU, U, ctx->UV, INSERT_VALUES, SCATTER_FORWARD);
135:   SNESSolve(ctx->snes, NULL, ctx->V);
136:   VecScatterBegin(ctx->scatterV, ctx->V, ctx->UV, INSERT_VALUES, SCATTER_FORWARD);
137:   VecScatterEnd(ctx->scatterV, ctx->V, ctx->UV, INSERT_VALUES, SCATTER_FORWARD);
138:   (*ctx->f)(t, ctx->UV, F);
139:   return 0;
140: }

142: /*
143:    Defines the nonlinear function that is passed to the nonlinear solver

145: */
146: PetscErrorCode SNESFunction(SNES snes, Vec V, Vec F, void *actx)
147: {
148:   AppCtx *ctx = (AppCtx *)actx;

151:   VecScatterBegin(ctx->scatterV, V, ctx->UV, INSERT_VALUES, SCATTER_FORWARD);
152:   VecScatterEnd(ctx->scatterV, V, ctx->UV, INSERT_VALUES, SCATTER_FORWARD);
153:   (*ctx->F)(ctx->t, ctx->UV, F);
154:   return 0;
155: }

157: /*TEST

159:     test:
160:       args: -ts_monitor -ts_view

162: TEST*/