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;

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

 85:   VecSet(UV,1.0);
 86:   TSSolve(ts,UV);
 87:   VecDestroy(&tsrhs);
 88:   VecDestroy(&UV);
 89:   TSDestroy(&ts);
 90:   PetscFinalize();
 91:   return 0;
 92: }

 94: /*
 95:    Defines the RHS function that is passed to the time-integrator.
 96: */
 97: PetscErrorCode TSFunctionRHS(TS ts,PetscReal t,Vec UV,Vec F,void *actx)
 98: {
 99:   AppCtx         *ctx = (AppCtx*)actx;

102:   VecSet(F,0.0);
103:   (*ctx->f)(t,UV,F);
104:   return 0;
105: }

107: /*
108:    Defines the nonlinear function that is passed to the time-integrator
109: */
110: PetscErrorCode TSFunctionI(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx)
111: {
112:   AppCtx         *ctx = (AppCtx*)actx;

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

120: /*TEST

122:     test:
123:       args:  -ts_view

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

129: TEST*/