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;

 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
 36: */
 37: PetscErrorCode F(PetscReal t, Vec UV, Vec F)
 38: {
 39:   const PetscScalar *u, *v;
 40:   PetscScalar       *f;
 41:   PetscInt           n, i;

 43:   PetscFunctionBeginUser;
 44:   PetscCall(VecGetLocalSize(UV, &n));
 45:   n = n / 2;
 46:   PetscCall(VecGetArrayRead(UV, &u));
 47:   v = u + n;
 48:   PetscCall(VecGetArrayWrite(F, &f));
 49:   for (i = 0; i < n; i++) f[i] = u[i] - v[i];
 50:   PetscCall(VecRestoreArrayRead(UV, &u));
 51:   PetscCall(VecRestoreArrayWrite(F, &f));
 52:   PetscFunctionReturn(PETSC_SUCCESS);
 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;

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

 89:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &ctx.snes));
 90:   PetscCall(SNESSetFromOptions(ctx.snes));
 91:   PetscCall(SNESSetFunction(ctx.snes, NULL, SNESFunction, &ctx));
 92:   PetscCall(SNESSetJacobian(ctx.snes, NULL, NULL, SNESComputeJacobianDefault, &ctx));
 93:   ctx.F = F;
 94:   PetscCall(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:   PetscCall(VecCreateMPI(PETSC_COMM_WORLD, 2, PETSC_DETERMINE, &ctx.UV));
 98:   i = 2 * rank;
 99:   PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 1, &i, PETSC_COPY_VALUES, &is));
100:   PetscCall(VecScatterCreate(U, NULL, ctx.UV, is, &ctx.scatterU));
101:   PetscCall(ISDestroy(&is));
102:   i = 2 * rank + 1;
103:   PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 1, &i, PETSC_COPY_VALUES, &is));
104:   PetscCall(VecScatterCreate(ctx.V, NULL, ctx.UV, is, &ctx.scatterV));
105:   PetscCall(ISDestroy(&is));

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

110:   PetscCall(VecDestroy(&ctx.V));
111:   PetscCall(VecDestroy(&ctx.UV));
112:   PetscCall(VecScatterDestroy(&ctx.scatterU));
113:   PetscCall(VecScatterDestroy(&ctx.scatterV));
114:   PetscCall(VecDestroy(&tsrhs));
115:   PetscCall(VecDestroy(&U));
116:   PetscCall(SNESDestroy(&ctx.snes));
117:   PetscCall(TSDestroy(&ts));
118:   PetscCall(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;

131:   PetscFunctionBeginUser;
132:   ctx->t = t;
133:   PetscCall(VecScatterBegin(ctx->scatterU, U, ctx->UV, INSERT_VALUES, SCATTER_FORWARD));
134:   PetscCall(VecScatterEnd(ctx->scatterU, U, ctx->UV, INSERT_VALUES, SCATTER_FORWARD));
135:   PetscCall(SNESSolve(ctx->snes, NULL, ctx->V));
136:   PetscCall(VecScatterBegin(ctx->scatterV, ctx->V, ctx->UV, INSERT_VALUES, SCATTER_FORWARD));
137:   PetscCall(VecScatterEnd(ctx->scatterV, ctx->V, ctx->UV, INSERT_VALUES, SCATTER_FORWARD));
138:   PetscCall((*ctx->f)(t, ctx->UV, F));
139:   PetscFunctionReturn(PETSC_SUCCESS);
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;

150:   PetscFunctionBeginUser;
151:   PetscCall(VecScatterBegin(ctx->scatterV, V, ctx->UV, INSERT_VALUES, SCATTER_FORWARD));
152:   PetscCall(VecScatterEnd(ctx->scatterV, V, ctx->UV, INSERT_VALUES, SCATTER_FORWARD));
153:   PetscCall((*ctx->F)(ctx->t, ctx->UV, F));
154:   PetscFunctionReturn(PETSC_SUCCESS);
155: }

157: /*TEST

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

162: TEST*/