Actual source code: ex3.c

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

\begin{eqnarray}
ys' = -2.0*\frac{-1.0+ys^2.0-\cos(t)}{2.0*ys}+0.05*\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf}-\frac{\sin(t)}{2.0*ys}\\
yf' = 0.05*\frac{-1.0+ys^2-\cos(t)}{2.0*ys}-\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf}-5.0*\frac{\sin(5.0*t)}{2.0*yf}\\
\end{eqnarray}

     12: #include <petscts.h>

     14: typedef struct {
     15:   PetscReal 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] = -2.0*(-1.0+u[0]*u[0]-PetscCosScalar(t))/(2.0*u[0])+0.05*(-2.0+u[1]*u[1]-PetscCosScalar(5.0*t))/(2.0*u[1])-PetscSinScalar(t)/(2.0*u[0]);
     28:   f[1] = 0.05*(-1.0+u[0]*u[0]-PetscCosScalar(t))/(2.0*u[0])-(-2.0+u[1]*u[1]-PetscCosScalar(5.0*t))/(2.0*u[1])-5.0*PetscSinScalar(5.0*t)/(2.0*u[1]);
     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] = -2.0*(-1.0+u[0]*u[0]-PetscCosScalar(t))/(2.0*u[0])+0.05*(-2.0+u[1]*u[1]-PetscCosScalar(5.0*t))/(2.0*u[1])-PetscSinScalar(t)/(2.0*u[0]);
     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] = 0.05*(-1.0+u[0]*u[0]-PetscCosScalar(t))/(2.0*u[0])-(-2.0+u[1]*u[1]-PetscCosScalar(5.0*t))/(2.0*u[1])-5.0*PetscSinScalar(5.0*t)/(2.0*u[1]);
     59:   VecRestoreArrayRead(U,&u);
     60:   VecRestoreArray(F,&f);
     61:   return(0);
     62: }

     64: static PetscErrorCode sol_true(PetscReal t,Vec U)
     65: {
     67:   PetscScalar    *u;

     70:   VecGetArray(U,&u);
     71:   u[0] = PetscSqrtScalar(1.0+PetscCosScalar(t));
     72:   u[1] = PetscSqrtScalar(2.0+PetscCosScalar(5.0*t));
     73:   VecRestoreArray(U,&u);
     74:   return(0);
     75: }

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

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

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

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

    119:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    120:     Set initial condition
    121:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    122:   VecGetArray(U,&u);
    123:   u[0] = PetscSqrtScalar(2.0);
    124:   u[1] = PetscSqrtScalar(3.0);
    125:   VecRestoreArray(U,&u);

    127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    128:      Create timestepping solver context
    129:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    130:   TSCreate(PETSC_COMM_WORLD,&ts);
    131:   TSSetType(ts,TSMPRK);

    133:   TSSetRHSFunction(ts,NULL,(TSRHSFunction)RHSFunction,&ctx);
    134:   TSRHSSplitSetIS(ts,"slow",iss);
    135:   TSRHSSplitSetIS(ts,"fast",isf);
    136:   TSRHSSplitSetRHSFunction(ts,"slow",NULL,(TSRHSFunction)RHSFunctionslow,&ctx);
    137:   TSRHSSplitSetRHSFunction(ts,"fast",NULL,(TSRHSFunction)RHSFunctionfast,&ctx);

    139:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    140:      Set initial conditions
    141:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    142:   TSSetSolution(ts,U);

    144:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    145:      Set solver options
    146:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    147:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ODE options","");
    148:   {
    149:     ctx.Tf = 0.3;
    150:     ctx.dt = 0.01;
    151:     PetscOptionsScalar("-Tf","","",ctx.Tf,&ctx.Tf,NULL);
    152:     PetscOptionsScalar("-dt","","",ctx.dt,&ctx.dt,NULL);
    153:   }
    154:   PetscOptionsEnd();
    155:   TSSetMaxTime(ts,ctx.Tf);
    156:   TSSetTimeStep(ts,ctx.dt);
    157:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
    158:   TSSetFromOptions(ts);

    160:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    161:      Solve linear system
    162:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    163:   TSSolve(ts,U);
    164:   VecView(U,PETSC_VIEWER_STDOUT_WORLD);

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

    174:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    175:      Print norm2 error
    176:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    177:   PetscPrintf(PETSC_COMM_WORLD,"l2 error norm: %g\n", error);

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

    193: /*TEST
    194:     build:
    195:       requires: !complex

    197:     test:

    199: TEST*/