Actual source code: ex9.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 ARKIMEX integrator on the entire DAE
 10: */

 12: /*
 13:    f(U,V) = U + V

 15: */
 16: PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F)
 17: {
 19:   VecWAXPY(F,1.0,U,V);
 20:   return 0;
 21: }

 23: /*
 24:    F(U,V) = U - V

 26: */
 27: PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F)
 28: {
 30:   VecWAXPY(F,-1.0,V,U);
 31:   return 0;
 32: }

 34: typedef struct {
 35:   Vec            U,V;
 36:   Vec            UF,VF;
 37:   VecScatter     scatterU,scatterV;
 38:   PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec);
 39:   PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec);
 40: } AppCtx;

 42: extern PetscErrorCode TSFunctionRHS(TS,PetscReal,Vec,Vec,void*);
 43: extern PetscErrorCode TSFunctionI(TS,PetscReal,Vec,Vec,Vec,void*);

 45: int main(int argc,char **argv)
 46: {
 47:   AppCtx         ctx;
 48:   TS             ts;
 49:   Vec            tsrhs,UV;
 50:   IS             is;
 51:   PetscInt       I;
 52:   PetscMPIInt    rank;

 54:   PetscInitialize(&argc,&argv,(char*)0,help);
 55:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 56:   TSCreate(PETSC_COMM_WORLD,&ts);
 57:   TSSetProblemType(ts,TS_NONLINEAR);
 58:   TSSetType(ts,TSROSW);
 59:   TSSetFromOptions(ts);
 60:   VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&tsrhs);
 61:   VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&UV);
 62:   TSSetRHSFunction(ts,tsrhs,TSFunctionRHS,&ctx);
 63:   TSSetIFunction(ts,NULL,TSFunctionI,&ctx);
 64:   ctx.f = f;
 65:   ctx.F = F;

 67:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.U);
 68:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.V);
 69:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.UF);
 70:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.VF);
 71:   I    = 2*rank;
 72:   ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is);
 73:   VecScatterCreate(ctx.U,NULL,UV,is,&ctx.scatterU);
 74:   ISDestroy(&is);
 75:   I    = 2*rank + 1;
 76:   ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is);
 77:   VecScatterCreate(ctx.V,NULL,UV,is,&ctx.scatterV);
 78:   ISDestroy(&is);

 80:   VecSet(UV,1.0);
 81:   TSSolve(ts,UV);
 82:   VecDestroy(&tsrhs);
 83:   VecDestroy(&UV);
 84:   VecDestroy(&ctx.U);
 85:   VecDestroy(&ctx.V);
 86:   VecDestroy(&ctx.UF);
 87:   VecDestroy(&ctx.VF);
 88:   VecScatterDestroy(&ctx.scatterU);
 89:   VecScatterDestroy(&ctx.scatterV);
 90:   TSDestroy(&ts);
 91:   PetscFinalize();
 92:   return 0;
 93: }

 95: /*
 96:    Defines the RHS function that is passed to the time-integrator.

 98: */
 99: PetscErrorCode TSFunctionRHS(TS ts,PetscReal t,Vec UV,Vec F,void *actx)
100: {
101:   AppCtx         *ctx = (AppCtx*)actx;

104:   VecSet(F,0.0);
105:   VecScatterBegin(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE);
106:   VecScatterEnd(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE);
107:   VecScatterBegin(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE);
108:   VecScatterEnd(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE);
109:   (*ctx->f)(t,ctx->U,ctx->V,ctx->UF);
110:   VecScatterBegin(ctx->scatterU,ctx->UF,F,INSERT_VALUES,SCATTER_FORWARD);
111:   VecScatterEnd(ctx->scatterU,ctx->UF,F,INSERT_VALUES,SCATTER_FORWARD);
112:   return 0;
113: }

115: /*
116:    Defines the nonlinear function that is passed to the time-integrator

118: */
119: PetscErrorCode TSFunctionI(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx)
120: {
121:   AppCtx         *ctx = (AppCtx*)actx;

124:   VecCopy(UVdot,F);
125:   VecScatterBegin(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE);
126:   VecScatterEnd(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE);
127:   VecScatterBegin(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE);
128:   VecScatterEnd(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE);
129:   (*ctx->F)(t,ctx->U,ctx->V,ctx->VF);
130:   VecScatterBegin(ctx->scatterV,ctx->VF,F,INSERT_VALUES,SCATTER_FORWARD);
131:   VecScatterEnd(ctx->scatterV,ctx->VF,F,INSERT_VALUES,SCATTER_FORWARD);
132:   return 0;
133: }