Actual source code: ex1.c

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

\begin{eqnarray}
ys' = \frac{ys}{a}\\
yf' = ys*cos(bt)\\
\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:   PetscErrorCode    ierr;
     21:   const PetscScalar *u;
     22:   PetscScalar       *f;

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

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

     41:   VecGetArrayRead(U,&u);
     42:   VecGetArray(F,&f);
     43:   f[0] = u[0]/ctx->a;
     44:   VecRestoreArrayRead(U,&u);
     45:   VecRestoreArray(F,&f);
     46:   return(0);
     47: }

     49: static PetscErrorCode RHSFunctionfast(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
     50: {
     51:   PetscErrorCode    ierr;
     52:   const PetscScalar *u;
     53:   PetscScalar       *f;

     56:   VecGetArrayRead(U,&u);
     57:   VecGetArray(F,&f);
     58:   f[0] = u[0]*PetscCosScalar(t*ctx->b);
     59:   VecRestoreArrayRead(U,&u);
     60:   VecRestoreArray(F,&f);
     61:   return(0);
     62: }

     64: /*
     65:   Define the analytic solution for check method easily
     66: */
     67: static PetscErrorCode sol_true(PetscReal t, Vec U,AppCtx *ctx)
     68: {
     70:   PetscScalar    *u;

     73:   VecGetArray(U,&u);
     74:   u[0] = PetscExpScalar(t/ctx->a);
     75:   u[1] = (ctx->a*PetscCosScalar(ctx->b*t)+ctx->a*ctx->a*ctx->b*PetscSinScalar(ctx->b*t))*PetscExpScalar(t/ctx->a)/(1+ctx->a*ctx->a*ctx->b*ctx->b);
     76:   VecRestoreArray(U,&u);
     77:   return(0);
     78: }

     80: int main(int argc,char **argv)
     81: {
     82:   TS             ts; /* ODE integrator */
     83:   Vec            U;  /* solution will be stored here */
     84:   Vec            Utrue;
     86:   PetscMPIInt    size;
     87:   AppCtx         ctx;
     88:   PetscScalar    *u;
     89:   IS             iss;
     90:   IS             isf;
     91:   PetscInt       *indicess;
     92:   PetscInt       *indicesf;
     93:   PetscInt       n=2;
     94:   PetscReal      error,tt;

     96:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     97:      Initialize program
     98:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     99:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
    100:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
    101:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");

    103:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    104:     Create index for slow part and fast part
    105:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    106:   PetscMalloc1(1,&indicess);
    107:   indicess[0] = 0;
    108:   PetscMalloc1(1,&indicesf);
    109:   indicesf[0] = 1;
    110:   ISCreateGeneral(PETSC_COMM_SELF,1,indicess,PETSC_COPY_VALUES,&iss);
    111:   ISCreateGeneral(PETSC_COMM_SELF,1,indicesf,PETSC_COPY_VALUES,&isf);

    113:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    114:     Create necesary vector
    115:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    116:   VecCreate(PETSC_COMM_WORLD,&U);
    117:   VecSetSizes(U,n,PETSC_DETERMINE);
    118:   VecSetFromOptions(U);
    119:   VecDuplicate(U,&Utrue);
    120:   VecCopy(U,Utrue);

    122:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    123:     Set runtime options
    124:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    125:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ODE options","");
    126:   {
    127:     ctx.a  = 2.0;
    128:     ctx.b  = 25.0;
    129:     PetscOptionsReal("-a","","",ctx.a,&ctx.a,NULL);
    130:     PetscOptionsReal("-b","","",ctx.b,&ctx.b,NULL);
    131:     ctx.Tf = 2;
    132:     ctx.dt = 0.01;
    133:     PetscOptionsReal("-Tf","","",ctx.Tf,&ctx.Tf,NULL);
    134:     PetscOptionsReal("-dt","","",ctx.dt,&ctx.dt,NULL);
    135:   }
    136:  PetscOptionsEnd();

    138:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    139:      Initialize the solution
    140:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    141:   VecGetArray(U,&u);
    142:   u[0] = 1.0;
    143:   u[1] = ctx.a/(1+ctx.a*ctx.a*ctx.b*ctx.b);
    144:   VecRestoreArray(U,&u);

    146:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    147:      Create timestepping solver context
    148:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    149:   TSCreate(PETSC_COMM_WORLD,&ts);
    150:   TSSetType(ts,TSMPRK);

    152:   TSSetRHSFunction(ts,NULL,(TSRHSFunction)RHSFunction,&ctx);
    153:   TSRHSSplitSetIS(ts,"slow",iss);
    154:   TSRHSSplitSetIS(ts,"fast",isf);
    155:   TSRHSSplitSetRHSFunction(ts,"slow",NULL,(TSRHSFunction)RHSFunctionslow,&ctx);
    156:   TSRHSSplitSetRHSFunction(ts,"fast",NULL,(TSRHSFunction)RHSFunctionfast,&ctx);

    158:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    159:      Set initial conditions
    160:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    161:   TSSetSolution(ts,U);

    163:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    164:      Set solver options
    165:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    166:   TSSetMaxTime(ts,ctx.Tf);
    167:   TSSetTimeStep(ts,ctx.dt);
    168:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
    169:   TSSetFromOptions(ts);

    171:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    172:      Solve linear system
    173:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    174:   TSSolve(ts,U);
    175:   VecView(U,PETSC_VIEWER_STDOUT_WORLD);

    177:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    178:      Check the error of the Petsc solution
    179:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    180:   TSGetTime(ts,&tt);
    181:   sol_true(tt,Utrue,&ctx);
    182:   VecAXPY(Utrue,-1.0,U);
    183:   VecNorm(Utrue,NORM_2,&error);

    185:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    186:      Print norm2 error
    187:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    188:   PetscPrintf(PETSC_COMM_WORLD,"l2 error norm = %g\n", error);

    190:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    191:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
    192:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    193:   VecDestroy(&U);
    194:   VecDestroy(&Utrue);
    195:   TSDestroy(&ts);
    196:   ISDestroy(&iss);
    197:   ISDestroy(&isf);
    198:   PetscFree(indicess);
    199:   PetscFree(indicesf);
    200:   PetscFinalize();
    201:   return ierr;
    202: }

    204: /*TEST
    205:     build:
    206:       requires: !complex

    208:     test:

    210: TEST*/