Actual source code: ex2.c

  1: static char help[] = "Basic problem for multi-rate method.\n";


\begin{eqnarray}
ys' = -sin(a*t)\\
yf' = bcos(b*t)ys-sin(b*t)sin(a*t)\\
\end{eqnarray}

 12: #include <petscts.h>

 14: typedef struct {
 15:   PetscReal a,b,Tf,dt;
 16: } AppCtx;

 18: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
 19: {
 20:   const PetscScalar *u;
 21:   PetscScalar       *f;

 23:   VecGetArrayRead(U,&u);
 24:   VecGetArray(F,&f);
 25:   f[0] = -PetscSinScalar(ctx->a*t);
 26:   f[1] = ctx->b*PetscCosScalar(ctx->b*t)*u[0]-PetscSinScalar(ctx->b*t)*PetscSinScalar(ctx->a*t);
 27:   VecRestoreArrayRead(U,&u);
 28:   VecRestoreArray(F,&f);
 29:   return 0;
 30:  }

 32: static PetscErrorCode RHSFunctionslow(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
 33: {
 34:   const PetscScalar *u;
 35:   PetscScalar       *f;

 37:   VecGetArrayRead(U,&u);
 38:   VecGetArray(F,&f);
 39:   f[0] = -PetscSinScalar(ctx->a*t);
 40:   VecRestoreArrayRead(U,&u);
 41:   VecRestoreArray(F,&f);
 42:   return 0;
 43: }

 45: static PetscErrorCode RHSFunctionfast(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
 46: {
 47:   const PetscScalar *u;
 48:   PetscScalar       *f;

 50:   VecGetArrayRead(U,&u);
 51:   VecGetArray(F,&f);
 52:   f[0] = ctx->b*PetscCosScalar(ctx->b*t)*u[0]-PetscSinScalar(ctx->b*t)*PetscSinScalar(ctx->a*t);
 53:   VecRestoreArrayRead(U,&u);
 54:   VecRestoreArray(F,&f);
 55:   return 0;
 56: }

 58: /*
 59:   Define the analytic solution for check method easily
 60: */
 61: static PetscErrorCode sol_true(PetscReal t,Vec U,AppCtx *ctx)
 62: {
 63:   PetscScalar    *u;

 65:   VecGetArray(U,&u);
 66:   u[0] = PetscCosScalar(ctx->a*t)/ctx->a;
 67:   u[1] = PetscSinScalar(ctx->b*t)*PetscCosScalar(ctx->a*t)/ctx->a;
 68:   VecRestoreArray(U,&u);
 69:   return 0;
 70: }

 72: int main(int argc,char **argv)
 73: {
 74:   TS             ts;            /* ODE integrator */
 75:   Vec            U;             /* solution will be stored here */
 76:   Vec            Utrue;
 78:   PetscMPIInt    size;
 79:   AppCtx         ctx;
 80:   PetscScalar    *u;
 81:   IS             iss;
 82:   IS             isf;
 83:   PetscInt       *indicess;
 84:   PetscInt       *indicesf;
 85:   PetscInt       n=2;
 86:   PetscScalar    error,tt;

 88:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 89:      Initialize program
 90:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 91:   PetscInitialize(&argc,&argv,(char*)0,help);
 92:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 95:  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 96:     Create index for slow part and fast part
 97:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 98:   PetscMalloc1(1,&indicess);
 99:   indicess[0]=0;
100:   PetscMalloc1(1,&indicesf);
101:   indicesf[0]=1;
102:   ISCreateGeneral(PETSC_COMM_SELF,1,indicess,PETSC_COPY_VALUES,&iss);
103:   ISCreateGeneral(PETSC_COMM_SELF,1,indicesf,PETSC_COPY_VALUES,&isf);

105:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106:     Create necessary vector
107:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108:   VecCreate(PETSC_COMM_WORLD,&U);
109:   VecSetSizes(U,n,PETSC_DETERMINE);
110:   VecSetFromOptions(U);
111:   VecDuplicate(U,&Utrue);
112:   VecCopy(U,Utrue);

114:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115:     Set runtime options
116:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ODE options","");
118:   {
119:     ctx.a  = 1.0;
120:     ctx.b  = 25.0;
121:     PetscOptionsScalar("-a","","",ctx.a,&ctx.a,NULL);
122:     PetscOptionsScalar("-b","","",ctx.b,&ctx.b,NULL);
123:     ctx.Tf = 5.0;
124:     ctx.dt = 0.01;
125:     PetscOptionsScalar("-Tf","","",ctx.Tf,&ctx.Tf,NULL);
126:     PetscOptionsScalar("-dt","","",ctx.dt,&ctx.dt,NULL);
127:   }
128:   PetscOptionsEnd();

130:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131:      Initialize the solution
132:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133:   VecGetArray(U,&u);
134:   u[0] = 1.0/ctx.a;
135:   u[1] = 0.0;
136:   VecRestoreArray(U,&u);

138:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139:      Create timestepping solver context
140:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141:   TSCreate(PETSC_COMM_WORLD,&ts);
142:   TSSetType(ts,TSMPRK);

144:   TSSetRHSFunction(ts,NULL,(TSRHSFunction)RHSFunction,&ctx);
145:   TSRHSSplitSetIS(ts,"slow",iss);
146:   TSRHSSplitSetIS(ts,"fast",isf);
147:   TSRHSSplitSetRHSFunction(ts,"slow",NULL,(TSRHSFunction)RHSFunctionslow,&ctx);
148:   TSRHSSplitSetRHSFunction(ts,"fast",NULL,(TSRHSFunction)RHSFunctionfast,&ctx);

150:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151:      Set initial conditions
152:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153:   TSSetSolution(ts,U);

155:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156:      Set solver options
157:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158:   TSSetMaxTime(ts,ctx.Tf);
159:   TSSetTimeStep(ts,ctx.dt);
160:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
161:   TSSetFromOptions(ts);

163:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164:      Solve linear system
165:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166:   TSSolve(ts,U);
167:   VecView(U,PETSC_VIEWER_STDOUT_WORLD);

169:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
170:      Check the error of the Petsc solution
171:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
172:   TSGetTime(ts,&tt);
173:   sol_true(tt,Utrue,&ctx);
174:   VecAXPY(Utrue,-1.0,U);
175:   VecNorm(Utrue,NORM_2,&error);

177:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178:      Print norm2 error
179:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180:   PetscPrintf(PETSC_COMM_WORLD,"l2 error norm: %g\n", error);

182:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
184:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185:   VecDestroy(&U);
186:   TSDestroy(&ts);
187:   VecDestroy(&Utrue);
188:   ISDestroy(&iss);
189:   ISDestroy(&isf);
190:   PetscFree(indicess);
191:   PetscFree(indicesf);
192:   PetscFinalize();
193:   return 0;
194: }

196: /*TEST
197:     build:
198:       requires: !complex

200:     test:

202: TEST*/