Actual source code: ex1.c

petsc-3.13.4 2020-08-01
Report Typos and Errors
  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 c99

208:     test:

210: TEST*/