Actual source code: ex2.c

  1: /*
  2:        Formatted test for TS routines.

  4:           Solves U_t=F(t,u)
  5:           Where:

  7:                   [2*u1+u2
  8:           F(t,u)= [u1+2*u2+u3
  9:                   [   u2+2*u3
 10:        We can compare the solutions from euler, beuler and SUNDIALS to
 11:        see what is the difference.

 13: */

 15: static char help[] = "Solves a linear ODE. \n\n";

 17: #include <petscts.h>
 18: #include <petscpc.h>

 20: extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
 21: extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);
 22: extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);
 23: extern PetscErrorCode Initial(Vec,void*);
 24: extern PetscErrorCode MyMatMult(Mat,Vec,Vec);

 26: extern PetscReal solx(PetscReal);
 27: extern PetscReal soly(PetscReal);
 28: extern PetscReal solz(PetscReal);

 30: int main(int argc,char **argv)
 31: {
 32:   PetscInt       time_steps = 100,steps;
 33:   Vec            global;
 34:   PetscReal      dt,ftime;
 35:   TS             ts;
 36:   Mat            A = 0,S;

 38:   PetscInitialize(&argc,&argv,(char*)0,help);
 39:   PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL);

 41:   /* set initial conditions */
 42:   VecCreate(PETSC_COMM_WORLD,&global);
 43:   VecSetSizes(global,PETSC_DECIDE,3);
 44:   VecSetFromOptions(global);
 45:   Initial(global,NULL);

 47:   /* make timestep context */
 48:   TSCreate(PETSC_COMM_WORLD,&ts);
 49:   TSSetProblemType(ts,TS_NONLINEAR);
 50:   TSMonitorSet(ts,Monitor,NULL,NULL);
 51:   dt = 0.001;

 53:   /*
 54:     The user provides the RHS and Jacobian
 55:   */
 56:   TSSetRHSFunction(ts,NULL,RHSFunction,NULL);
 57:   MatCreate(PETSC_COMM_WORLD,&A);
 58:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3);
 59:   MatSetFromOptions(A);
 60:   MatSetUp(A);
 61:   RHSJacobian(ts,0.0,global,A,A,NULL);
 62:   TSSetRHSJacobian(ts,A,A,RHSJacobian,NULL);

 64:   MatCreateShell(PETSC_COMM_WORLD,3,3,3,3,NULL,&S);
 65:   MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MyMatMult);
 66:   TSSetRHSJacobian(ts,S,A,RHSJacobian,NULL);

 68:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
 69:   TSSetFromOptions(ts);

 71:   TSSetTimeStep(ts,dt);
 72:   TSSetMaxSteps(ts,time_steps);
 73:   TSSetMaxTime(ts,1);
 74:   TSSetSolution(ts,global);

 76:   TSSolve(ts,global);
 77:   TSGetSolveTime(ts,&ftime);
 78:   TSGetStepNumber(ts,&steps);

 80:   /* free the memories */

 82:   TSDestroy(&ts);
 83:   VecDestroy(&global);
 84:   MatDestroy(&A);
 85:   MatDestroy(&S);

 87:   PetscFinalize();
 88:   return 0;
 89: }

 91: PetscErrorCode MyMatMult(Mat S,Vec x,Vec y)
 92: {
 93:   const PetscScalar  *inptr;
 94:   PetscScalar        *outptr;

 96:   VecGetArrayRead(x,&inptr);
 97:   VecGetArrayWrite(y,&outptr);

 99:   outptr[0] = 2.0*inptr[0]+inptr[1];
100:   outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2];
101:   outptr[2] = inptr[1]+2.0*inptr[2];

103:   VecRestoreArrayRead(x,&inptr);
104:   VecRestoreArrayWrite(y,&outptr);
105:   return 0;
106: }

108: /* this test problem has initial values (1,1,1).                      */
109: PetscErrorCode Initial(Vec global,void *ctx)
110: {
111:   PetscScalar    *localptr;
112:   PetscInt       i,mybase,myend,locsize;

114:   /* determine starting point of each processor */
115:   VecGetOwnershipRange(global,&mybase,&myend);
116:   VecGetLocalSize(global,&locsize);

118:   /* Initialize the array */
119:   VecGetArrayWrite(global,&localptr);
120:   for (i=0; i<locsize; i++) localptr[i] = 1.0;

122:   if (mybase == 0) localptr[0]=1.0;

124:   VecRestoreArrayWrite(global,&localptr);
125:   return 0;
126: }

128: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
129: {
130:   VecScatter        scatter;
131:   IS                from,to;
132:   PetscInt          i,n,*idx;
133:   Vec               tmp_vec;
134:   PetscErrorCode    ierr;
135:   const PetscScalar *tmp;

137:   /* Get the size of the vector */
138:   VecGetSize(global,&n);

140:   /* Set the index sets */
141:   PetscMalloc1(n,&idx);
142:   for (i=0; i<n; i++) idx[i]=i;

144:   /* Create local sequential vectors */
145:   VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);

147:   /* Create scatter context */
148:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
149:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
150:   VecScatterCreate(global,from,tmp_vec,to,&scatter);
151:   VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
152:   VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);

154:   VecGetArrayRead(tmp_vec,&tmp);
155:   PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e u = %14.6e  %14.6e  %14.6e \n",
156:                      (double)time,(double)PetscRealPart(tmp[0]),(double)PetscRealPart(tmp[1]),(double)PetscRealPart(tmp[2]));
157:   PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e errors = %14.6e  %14.6e  %14.6e \n",
158:                      (double)time,(double)PetscRealPart(tmp[0]-solx(time)),(double)PetscRealPart(tmp[1]-soly(time)),(double)PetscRealPart(tmp[2]-solz(time)));
159:   VecRestoreArrayRead(tmp_vec,&tmp);
160:   VecScatterDestroy(&scatter);
161:   ISDestroy(&from);
162:   ISDestroy(&to);
163:   PetscFree(idx);
164:   VecDestroy(&tmp_vec);
165:   return 0;
166: }

168: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
169: {
170:   PetscScalar       *outptr;
171:   const PetscScalar *inptr;
172:   PetscInt          i,n,*idx;
173:   IS                from,to;
174:   VecScatter        scatter;
175:   Vec               tmp_in,tmp_out;

177:   /* Get the length of parallel vector */
178:   VecGetSize(globalin,&n);

180:   /* Set the index sets */
181:   PetscMalloc1(n,&idx);
182:   for (i=0; i<n; i++) idx[i]=i;

184:   /* Create local sequential vectors */
185:   VecCreateSeq(PETSC_COMM_SELF,n,&tmp_in);
186:   VecDuplicate(tmp_in,&tmp_out);

188:   /* Create scatter context */
189:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
190:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
191:   VecScatterCreate(globalin,from,tmp_in,to,&scatter);
192:   VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
193:   VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
194:   VecScatterDestroy(&scatter);

196:   /*Extract income array */
197:   VecGetArrayRead(tmp_in,&inptr);

199:   /* Extract outcome array*/
200:   VecGetArrayWrite(tmp_out,&outptr);

202:   outptr[0] = 2.0*inptr[0]+inptr[1];
203:   outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2];
204:   outptr[2] = inptr[1]+2.0*inptr[2];

206:   VecRestoreArrayRead(tmp_in,&inptr);
207:   VecRestoreArrayWrite(tmp_out,&outptr);

209:   VecScatterCreate(tmp_out,from,globalout,to,&scatter);
210:   VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
211:   VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);

213:   /* Destroy idx aand scatter */
214:   ISDestroy(&from);
215:   ISDestroy(&to);
216:   VecScatterDestroy(&scatter);
217:   VecDestroy(&tmp_in);
218:   VecDestroy(&tmp_out);
219:   PetscFree(idx);
220:   return 0;
221: }

223: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat BB,void *ctx)
224: {
225:   PetscScalar       v[3];
226:   const PetscScalar *tmp;
227:   PetscInt          idx[3],i;

229:   idx[0]=0; idx[1]=1; idx[2]=2;
230:   VecGetArrayRead(x,&tmp);

232:   i    = 0;
233:   v[0] = 2.0; v[1] = 1.0; v[2] = 0.0;
234:   MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);

236:   i    = 1;
237:   v[0] = 1.0; v[1] = 2.0; v[2] = 1.0;
238:   MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);

240:   i    = 2;
241:   v[0] = 0.0; v[1] = 1.0; v[2] = 2.0;
242:   MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);

244:   MatAssemblyBegin(BB,MAT_FINAL_ASSEMBLY);
245:   MatAssemblyEnd(BB,MAT_FINAL_ASSEMBLY);

247:   if (A != BB) {
248:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
249:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
250:   }
251:   VecRestoreArrayRead(x,&tmp);

253:   return 0;
254: }

256: /*
257:       The exact solutions
258: */
259: PetscReal solx(PetscReal t)
260: {
261:   return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0)) +
262:          PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0));
263: }

265: PetscReal soly(PetscReal t)
266: {
267:   return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/PetscSqrtReal(2.0) +
268:          PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/PetscSqrtReal(2.0);
269: }

271: PetscReal solz(PetscReal t)
272: {
273:   return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0)) +
274:          PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0));
275: }

277: /*TEST

279:     test:
280:       suffix: euler
281:       args: -ts_type euler
282:       requires: !single

284:     test:
285:       suffix: beuler
286:       args:   -ts_type beuler
287:       requires: !single

289: TEST*/