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: */


 15: /*
 16:    f(U,V) = U + V

 18: */
 19: PetscErrorCode f(PetscReal t,Vec UV,Vec F)
 20: {
 21:   PetscErrorCode    ierr;
 22:   const PetscScalar *u,*v;
 23:   PetscScalar       *f;
 24:   PetscInt          n,i;

 27:   VecGetLocalSize(UV,&n);
 28:   n    = n/2;
 29:   VecGetArrayRead(UV,&u);
 30:   v    = u + n;
 31:   VecGetArray(F,&f);
 32:   for (i=0; i<n; i++) f[i] = u[i] + v[i];
 33:   VecRestoreArrayRead(UV,&u);
 34:   VecRestoreArray(F,&f);
 35:   return(0);
 36: }

 40: /*
 41:    F(U,V) = U - V

 43: */
 44: PetscErrorCode F(PetscReal t,Vec UV,Vec F)
 45: {
 46:   PetscErrorCode    ierr;
 47:   const PetscScalar *u,*v;
 48:   PetscScalar       *f;
 49:   PetscInt          n,i;

 52:   VecGetLocalSize(UV,&n);
 53:   n    = n/2;
 54:   VecGetArrayRead(UV,&u);
 55:   v    = u + n;
 56:   VecGetArray(F,&f);
 57:   for (i=0; i<n; i++) f[i] = u[i] - v[i];
 58:   VecRestoreArrayRead(UV,&u);
 59:   VecRestoreArray(F,&f);
 60:   return(0);
 61: }

 63: typedef struct {
 64:   PetscReal      t;
 65:   SNES           snes;
 66:   Vec            UV,V;
 67:   VecScatter     scatterU,scatterV;
 68:   PetscErrorCode (*f)(PetscReal,Vec,Vec);
 69:   PetscErrorCode (*F)(PetscReal,Vec,Vec);
 70: } AppCtx;

 72: extern PetscErrorCode TSFunction(TS,PetscReal,Vec,Vec,void*);
 73: extern PetscErrorCode SNESFunction(SNES,Vec,Vec,void*);

 77: int main(int argc,char **argv)
 78: {
 80:   AppCtx         ctx;
 81:   TS             ts;
 82:   Vec            tsrhs,U;
 83:   IS             is;
 84:   PetscInt       I;
 85:   PetscMPIInt    rank;

 87:   PetscInitialize(&argc,&argv,(char*)0,help);
 88:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 89:   TSCreate(PETSC_COMM_WORLD,&ts);
 90:   TSSetProblemType(ts,TS_NONLINEAR);
 91:   TSSetType(ts,TSEULER);
 92:   TSSetFromOptions(ts);
 93:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&tsrhs);
 94:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&U);
 95:   TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx);
 96:   ctx.f = f;

 98:   SNESCreate(PETSC_COMM_WORLD,&ctx.snes);
 99:   SNESSetFromOptions(ctx.snes);
100:   SNESSetFunction(ctx.snes,NULL,SNESFunction,&ctx);
101:   SNESSetJacobian(ctx.snes,NULL,NULL,SNESComputeJacobianDefault,&ctx);
102:   ctx.F = F;
103:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.V);

105:   /* Create scatters to move between separate U and V representation and UV representation of solution */
106:   VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&ctx.UV);
107:   I    = 2*rank;
108:   ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is);
109:   VecScatterCreate(U,NULL,ctx.UV,is,&ctx.scatterU);
110:   ISDestroy(&is);
111:   I    = 2*rank + 1;
112:   ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is);
113:   VecScatterCreate(ctx.V,NULL,ctx.UV,is,&ctx.scatterV);
114:   ISDestroy(&is);

116:   VecSet(U,1.0);
117:   TSSolve(ts,U);

119:   VecDestroy(&ctx.V);
120:   VecDestroy(&ctx.UV);
121:   VecScatterDestroy(&ctx.scatterU);
122:   VecScatterDestroy(&ctx.scatterV);
123:   VecDestroy(&tsrhs);
124:   VecDestroy(&U);
125:   SNESDestroy(&ctx.snes);
126:   TSDestroy(&ts);
127:   PetscFinalize();
128:   return 0;
129: }

133: /*
134:    Defines the RHS function that is passed to the time-integrator. 

136:    Solves F(U,V) for V and then computes f(U,V)

138: */
139: PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
140: {
141:   AppCtx         *ctx = (AppCtx*)actx;

145:   ctx->t = t;
146:   VecScatterBegin(ctx->scatterU,U,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
147:   VecScatterEnd(ctx->scatterU,U,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
148:   SNESSolve(ctx->snes,NULL,ctx->V);
149:   VecScatterBegin(ctx->scatterV,ctx->V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
150:   VecScatterEnd(ctx->scatterV,ctx->V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
151:   (*ctx->f)(t,ctx->UV,F);
152:   return(0);
153: }

157: /*
158:    Defines the nonlinear function that is passed to the nonlinear solver

160: */
161: PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx)
162: {
163:   AppCtx         *ctx = (AppCtx*)actx;

167:   VecScatterBegin(ctx->scatterV,V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
168:   VecScatterEnd(ctx->scatterV,V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
169:   (*ctx->F)(ctx->t,ctx->UV,F);
170:   return(0);
171: }