Actual source code: ex3.c

petsc-3.12.0 2019-09-29
Report Typos and Errors
  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,(TSRHSFunctionslow)RHSFunctionslow,&ctx);
137:   TSRHSSplitSetRHSFunction(ts,"fast",NULL,(TSRHSFunctionfast)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: }