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:   PetscErrorCode    ierr;
 19:   const PetscScalar *u,*v;
 20:   PetscScalar       *f;
 21:   PetscInt          n,i;

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

 35: /*
 36:    F(U,V) = U - V
 37: */
 38: PetscErrorCode F(PetscReal t,Vec UV,Vec F)
 39: {
 40:   PetscErrorCode    ierr;
 41:   const PetscScalar *u,*v;
 42:   PetscScalar       *f;
 43:   PetscInt          n,i;

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

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

 66: extern PetscErrorCode TSFunction(TS,PetscReal,Vec,Vec,void*);
 67: extern PetscErrorCode SNESFunction(SNES,Vec,Vec,void*);

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

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

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

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

109:   VecSet(U,1.0);
110:   TSSolve(ts,U);

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

124: /*
125:    Defines the RHS function that is passed to the time-integrator.

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

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

145: /*
146:    Defines the nonlinear function that is passed to the nonlinear solver

148: */
149: PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx)
150: {
151:   AppCtx         *ctx = (AppCtx*)actx;

155:   VecScatterBegin(ctx->scatterV,V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
156:   VecScatterEnd(ctx->scatterV,V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
157:   (*ctx->F)(ctx->t,ctx->UV,F);
158:   return(0);
159: }

161: /*TEST

163:     test:
164:       args: -ts_monitor -ts_view

166: TEST*/